ON THE NUMERICAL EVALUATION OF DISTRIBUTIONS 
IN RANDOM MATRIX THEORY: A REVIEW 



FOLKMAR BORNEMANN 

Abstract. In this paper we review and compare the numerical evaluation of those 
probability distributions in random matrix theory that are analytically represented 
in terms of Painleve transcendents or Fredholm determinants. Concrete examples 
for the Gaussian and Laguerre (Wishart) /3-ensembles and their various scaling 
limits are discussed. We argue that the numerical approximation of Fredholm 
determinants is the conceptually more simple and efficient of the two approaches, 
I easily generalized to the computation of joint probabilities and correlations. Having 

^/-^ the means for extensive numerical explorations at hand, we discovered new and 

surprising determinantal formulae for the /c-th largest (or smallest) level in the 

I I edge scaling limits of the Orthogonal and Symplectic Ensembles; formulae that in 

turn led to improved numerical evaluations. The paper comes with a toolbox of 
Matlab functions that facilitates further mathematical experiments by the reader. 

-I— > 

C3 

i. Introduction 

Random Matrix Theory (RMT) has found many applications, most notably in 
physics, multivariate statistics, electrical engineering, and finance. As soon as there 
—i is the need for specific numbers, such as moments, quantiles, or correlations, the 

actual numerical evaluation of the underlying probability distributions becomes of 
interest. Without additional structure there would be, in general, only one method: 
Monte Carlo simulation. However, because of the universality of certain scaling 
limits (for a review see, e.g., Deift 2007), a family of distinguished distribution 
functions enters which is derived from highly structured matrix models enjoying 
closed analytic solutions. These functions constitute a new class of special functions 
comparable in import to the classic distributions of probability theory. This paper 
addresses the accurate numerical evaluation 1 of many of these functions on the 
one hand and shows, on the other hand, that such work facilitates numerical 
explorations that may lead, in the sense of Experimental Mathematics (Borwein 
and Bailey 2004), to new theoretical discoveries, see the results of Section 6. 

1.1. The Common Point of View. The closed analytic solutions alluded to above 
are based (for deeper reasons or because of contingency) on two concurrent tools: 
Fredholm determinants of integral operators and Painleve transcendents. Concern- 
ing the question which of them is better suited to be attacked numerically, there 
has been a prevailing point of view for the last 15 years or so, explicitly formulated 
by Tracy and Widom (2000, Footnote 10): "Without the Painleve representation, 



> 

X 



2000 Mathematics Subject Classification. Primary 15A52, 65R20; Secondary 33E17, 47G10. 
'Limiting the means, as customary in numerical analysis for reasons of efficiency and strict adher- 
ence to numerical stability, to IEEE double precision hardware arithmetic (about 16 digits precision). 



1 



2 



FOLKMAR BORNEMANN 



the numerical evaluation of the Fredholm determinants is quite involved." To 
understand the possible genesis of this point of view let us recall the results for 
the two most important scaling limits of the Gaussian Unitary Ensemble (GUE). 

1.1.1. Level Spacing Function of GUE. The large matrix limit of GUE, scaled for level 
spacing 1 in the bulk, yields the function 

E 2 (0;s) = P(no levels lie in (0,s)). (1.1) 

Gaudin (1961) showed that this function can be represented as a Fredholm deter- 
minant, namely, 2 

E 2 (0;s) = det - K sin \ L2{0 ^ , K sin (x,y) = sinc(7i(x - y)). (1.2) 

He proceeded by showing that the eigenfunctions of this selfadjoint integral 
operator are the radial prolate spheroidal wave functions with certain parameters. 
Using tables (Stratton, Morse, Chu, Little and Corbato 1956) of these special 
functions he was finally able to evaluate E 2 (0;s) numerically. 3 On the other hand, 
in an admirably intricate analytic tour de force Jimbo, Miwa, Mori and Sato (1980) 
expressed the Fredholm determinant by 

£ s (0;s) = exp (- j™ ^ dx\ (1.3) 

in terms of the Jimbo-Miwa-Okamoto cr-form of Painleve V, namely 

x x 2 

{xcr xx ) 2 = 4(<7 - xcr x )(xa x - a - 0%), < 7 ( x )~— + ^2 ( x ~* °)- 

1.1.2. The Tracy-Widom distribution. The large matrix limit of GUE, scaled for the 
fluctuations at the soft edge (that is, the maximum eigenvalue), yields the function 

F 2 (s) = P(no levels lie in (s,oo)). (1.5) 

Implicitly known for quite some time (see, e.g., Bronk 1964, Brezin and Kazakov 
1990, Moore 1990), the determinantal representation 

, . ^ f . \ „ , s Ai(x)Ai'(i/) - Ai'fx)Ai(v) , , 
F 2 (s) = det (I- K M \ L 2 (Sr00) ), K Ai (x,y) = ^^ V ^_ y K ' U - , (1.6) 

was spelt out by Forrester (1993) and by Tracy and Widom (1993b). The search 
for an analogue to Gaudin's method remained unsuccessful since there is no 
solution of the corresponding eigenvalue problem known in terms of classic 



2 We use the same symbol K to denote both, the integral operator K \x acting on the Hilbert space X 
and its kernel function K(x,y). 

^Strictly speaking Gaudin (1961) was concerned with evaluating the level spacing function Ei (0; s) 
of GOE that he represented as the Fredholm determinant of the even sine kernel, see (5.8a) below. 
However, the extension of Gaudin's method to £2(0; s) is fairly straightforward, see Dyson (1962) 
for a determinantal formula that is equivalent to (1.2), namely (5.7) with k = 0, and Kahn (1963) for 
subsequent numerical work. As was pointed out by Odlyzko (1987, p. 305), who himself had calculated 
E2(0;s) by Gaudin's method using Van Buren's implementation of the prolate wave functions, Kahn's 
tables, reproduced in the first 1967 edition of Mehta's book, are rather inaccurate. In contrast, the tables 
in Mehta and des Cloizeaux (1972), reproduced in the second 1991 and third 2002 edition of Mehta's 
book, are basically accurate with the proviso that the arguments have to be read not as the displayed 
four digit numbers but rather as s = 2t/n with t = 0.0, 0.1, 0.2, 0.3, etc. A modern implementation of 
Gaudin's method for £2(0; s), using Mathematica's fairly recent ability to evaluate the prolate wave 
functions, can be found in Bornemann (2010a, §7.1). 
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special functions (Mehta 2004, p. 453). It was therefore a major breakthrough when 
Tracy and Widom (1993b, 1994a) derived their now famous representation 

F2(s) = exp ^— J (x — s)u(x) 2 dx^j (1.7) 

in terms of the Hastings-McLeod (1980) solution u(x) of Painleve II, namely 

u xx = 2m 3 + xu, u{x) ~ Ai{x) (x — > 00). (1.8) 

Subsequent numerical evaluations were then, until the recent work of Bornemann 
(2010a), exclusively based on solving this asymptotic initial value problem. 



1.2. Challenging the Common Point of View. In this paper we challenge the 
common point of view that a Painleve representation would be, at least numerically, 
preferable to a Fredholm determinant formula. We do so from the following angles: 
Simplicity, efficiency, accuracy, and extendibility. Let us briefly indicate the rationale 
for our point of view. 

(1) The numerical evaluation of Painleve transcendents encountered in RMT 
is more involved as one would think at first sight. For reasons of numerical 
stability one needs additional, deep analytic knowledge, namely, asymp- 
totic expansions of the corresponding connection formulae (see Section 3). 

(2) There is an extremely simple, fast, accurate, and general numerical method 
for evaluating Fredholm determinants (see Section 4). 

(3) Multivariate functions such as joint probability distributions often have a 
representation by a Fredholm determinant (see Section 8). On the other 
hand, if available at all, a representation in terms of a nonlinear partial 
differential equation is of very limited numerical use right now. 



1.3. Outline of the Paper. In Section 2 we collect some fundamental functions 
of RMT whose numerical solutions will play a role in the sequel. The intricate 
issues of a numerical solution of the Painleve transcendents encountered in RMT 
are subject of Section 3. An exposition of Bornemann's (2010a) method for the 
numerical evaluation of Fredholm determinants is given in Section 4. The nu- 
merical evaluation of the A:-level spacings in the bulk of GOE, GUE, and GSE, by 
using Fredholm determinants, is addressed in Section 5. In Section 6 we get to 
new determinantal formulae for the distributions of the A:-th largest level in the 
soft edge scaling limit of GOE and GSE. These formulae rely on a determinantal 
identity that we found by extensive numerical experiments before proving it. By 
a powerful structural analogy, in Section 7 these formulae are easily extended to 
the A:-th smallest level in the hard edge scaling limit of LOE and LSE. In Section 8, 
we discuss some examples of joint probabilities, like the one for the largest two 
eigenvalues of GUE at the soft edge or the one of the Airy process for two different 
times. Finally, in Section 9, we give a short introduction into using the Matlab 
toolbox that comes with this paper. 
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2. Some Distribution Functions of RMT and Their Representation 

In this section we collect some fundamental functions of RMT whose numerical 
evaluation will explicitly be discussed in the sequel. We do not strive for com- 
pleteness here; but we will give sufficiently many examples to be able to judge of 
simplicity and generality of the numerical approaches later on. 

We confine ourselves to the Gaussian (Hermite) and Laguerre (Wishart) ensem- 
bles. That is, we consider n x n random matrix ensembles with real (nonnegative 
in the case of the Laguerre ensemble) spectrum such that the joint probability 
distribution of its (unordered) eigenvalues is given by 

p(xi,...,X n ) = C Y\v}{Xi) Y[\ X i - X jf- C 2 - 1 ) 

i i<j 

Here, /5 takes the values 1, 2, or 4 (Dyson's "three fold way"). The weight function 
w (x) will be either a Gaussian or, in the case of the Laguerre ensembles, a function 
of the type w a (x) = x a e~ x (a > —1). 

2.1. Gaussian Ensembles. Here, we take the Gaussian weight functions 4 

• w(x) = e~* 2 / 2 for 8 = 1, the Gaussian Orthogonal Ensemble (GOE), 

• w(x) = e~ x for f> — 2, the Gaussian Unitary Ensemble (GUE), 

2 - 

• wix) = e~ x for B = 4, the Gaussian Symplectic Ensemble (GSE). 
We define, for an open interval / C R, the basic quantity 

= P (exactly k eigenvalues of the n x n Gaussian S-ensemble lie in /) . (2.2) 
More general quantities will be considered in Section 8.3. 

2.1.1. Scaling limits. The bulk scaling limit is given by (see Mehta 2004, §§6.3, 7.2, 
and 11.7) 



Ef ulk \k;J)={" ■ (2.3) 



lim E { " } (k; n 2- 1/2 n- 1/2 }) P> = \oxP> = 2, 

lim E^\k;nn- y2 J) 8 = 4. 

n— 't<x> P 

The soft edge scaling limit is given by (see Forrester and Rains 2001, p. 194) 

f lim E { " ] (Jfc; V2n~ + 2- 1/2 n~ 1/6 J) B = 1 or B = 2, 

Ei oit) (k;J)= P (2.4) 

' \lunE^ 2 \k;V2n + 2- y2 n-^J) 6 = 4. 



4We follow Forrester and Rains (2001) in the choice of the variances of the Gaussian weights. Note 
that Mehta (2004, Chap. 3), such as Tracy and Widom in most of their work, uses w(x) = e - ^ 1 ' 2 ^ 2 . 
However, one has to be alert: from p. 175 onwards, Mehta uses w(x) = e for f> = 4 in his book, too. 
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2.1.2. Level spacing. The A:-level spacing function Ep(k;s) of Mehta (2004, §6.1.2) is 

E^(fc;s)=£< j bulk) (fc;(0 / s)). (2.5) 

The A:-level spacing density pp(k;s), that is, the probability density of the distance 
of a level in the bulk to its (k + l)-st next neighbor is (see Mehta 2004, Eq. (6.1.18)) 

Pp (k;s) = ^^(k + l-i)Ef i (r,s) (* = 0,1,2,...). (2.6) 

flS 7=0 

Since the bulk scaling limit was made such that the expected distance between 
neighboring eigenvalues is one, we have (see Mehta 2004, Eq. (6.1.26)) 



^ k;s) = l, spp(k;s)ds = fc + 1. (2.7) 

Likewise, there holds for s ^ 

00 00 

£E /5 (fc;s)=l / J^kEpfcs) =s; (2.8) 

fc=0 k=0 

see, e.g., Mehta (2004, Eqs. (6.1.19/20)) or Deift (1999b, p. 119). Such constraints 
are convenient means to assess the accuracy of numerical methods (see Exam- 
ples 9.3.2/9.4.2 and, for a generalization, Example 9.3.3). 

2.1.3. Distribution of the k-th largest eigenvalue. The cumulative distribution function 
of the A:-th largest eigenvalue is, in the soft edge scaling limit, 

^(^) = EEi soft) (j;(^)). (2.9) 

;=0 

The famous Tracy-Widom (1996) distributions F j g(s) are given by 

(F«(l;s) /3 = lor/3 = 2, 

F/3(s)= n (2.10) 

|^(1; \/2s) 0=4. 

2.2. Determinantal Representations for Gaussian Ensembles. 

2.2.1. GLfE. Here, the basic formula is (see Deift 1999b, §5.4) 



(2.11a) 



z=l 



dJ°(z;/) = det (1 -zK n \ L 2 i;) ) , (2.11b) 
with the Hermite kernel (the second form follows from Christoffel-Darboux) 

K n{x , y) = X>(^(y) = /| ^W»»-i(y)-» B -iW^(y) (2 . 12) 

that is built from the L 2 (R)-orthonormal system of the Hermite functions 

<p k (x) = *\ . (2.13) 
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The bulk scaling limit is given by (see Mehta 2004, §A.io) 



£r k) (^;/) = ^^Dr k) ( Z ;/) 



z=l 



D 



(bulk), . 



z;J) = det (1-zK 



v sinlL 2 (/) 



(2.14a) 



(2.14b) 



with the sine kernel K s ; n defined in (1.2). The soft edge scaling limit is given by 
(see Forrester 1993, §3.1) 



.(soft) 



(*;J) 



(-i) k d k 

M dz*" 2 



D, 



(soft) 



z;7) 



z=l 



D( soft) ( Z ;/)=det(l-zX Ai r L2(/) 
with the Airy kernel K^i defined in (1.6). 



(2.15a) 
(2.15b) 



2.2.2. GSE. Using Dyson's quaternion determinants one can find a determinantal 

formula for E^"' (k; /) which involves a finite rank matrix kernel (see Mehta 2004, 
Chap. 8) — a formula that is amenable to the numerical methods of Section 4. We 
confine ourselves to the soft edge scaling limit of this formula which yields (Tracy 
and Widom 2005) 



where the entries of the matrix kernel determinant 

s sd\ 



D 4 (z;/) 



z=l 



D 4 (z;J) =det I 



2 yis s* J 



\lHJ)(BLZ(J) 



(2.16a) 



(2.16b) 



are given by (with the adjoint kernel S* (x, y) = S{y,x) obtained from transposition) 



1 f°° 
S(x,y) = K Ai (x,y) - -Ai(x) J Ai(rj) dt], 



SD{x,y) = -d y K Ai (x,y) - -Ai(x)Ai(y), 



(2.17a) 
(2.17b) 



K Ai (5,y)dZ + -j x Ai(?) jf Ai(r,)d V . (2.17c) 

Albeit this expression is also amenable to the numerical methods of Section 4, we 
will discuss, for the significant special cases 



Ef ulk) (fc;(0,s)) and e£° ft) (fc; ( S ,oo)), 



(2.18) 



alternative determinantal formulae that are far more efficient numerically see 
Section 5 and Section 6, respectively. 
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2.2.3. GOE. There are also determinantal formulae for £j bulk ) j) anc [ its various 
scaling limits, see Mehta (2004, Chap. 7) and Tracy and Widom (2005). However, 
these determinantal formulae are based on matrix kernels that involve a discon- 
tinuous term. To be specific we recall the result for the soft edge scaling limit: 



where the entries of the matrix kernel determinant 



(2.19a) 

z=l 



Di(z;/) = det (l - z M \ Xl (j)®x 2 (j)j (2.19b) 

are given by (with the adjoint kernel S* (x, y) = S{y,x) obtained from transposition) 

S(x,y) = K Ai (x,y) + ^Ai(x) ^1 - J Ai(jy) dr^j , (2.20a) 

1 

SD(x,y) = -dyK Ai (x,y) - -Ai(x)Ai(y), (2.20b) 

/■CO 

IS £ {x,y) = -e(x - y) - / K M (£,y) dl, 

J X 

+ l - Ai(£) dl + ^ Ai(£) ^ Ai(jy) ^ . (2.20c) 
with the discontinuous function 

e(x) = 2sign(x). (2.2od) 

This discontinuity poses considerable difficulties for the proper theoretical justifi- 
cation of the operator determinant: for appropriately chosen weighted L 2 spaces 
Xj(J) and X2Q), the matrix kernel operator is a Hilbert-Schmidt operator with 
a trace class diagonal and the determinant has to be understood as a Hilbert- 
Carleman regularized determinant (see Tracy and Widom 2005, p. 2199). Moreover, 
it renders the unmodified numerical methods of Section 4 rather inefficient. Nev- 
ertheless, for the significant special cases 

£{ bulk) (fc;(0,s)) and E[ soit) {k; (s,co)) (2.21) 

there are alternative determinantal formulae, which are amenable to an efficient 
numerical evaluation, see Section 5 and Section 6, respectively. 

2.3. Painleve Representations for Gaussian Ensembles. For the important fam- 
ily of integrable kernels (see Deift 1999a), Tracy and Widom (1993a, 1993b, 1994a, 
1994c) found a general method to represent determinants of the form 



det [l- Z K\ L 2 { 

a b) ) (2.22) 

explicitly by a system of partial differential equations with respect to the indepen- 
dent variables a and b. Fixing one of the bounds yields an ordinary differential 
equation (that notwithstanding depends on the fixed bound). In RMT, this or- 
dinary differential equation turned out, case by case, to be a Painleve equation. 
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The typical choices of intervals with a fixed bound are / = (0, s) or / = (s, oo), 
depending on whether one looks at the bulk or the soft edge of the spectrum. 

2.3.1. GUE. Tracy and Widom (1994c) calculated, for the determinant (2.11b) with 
/ = (s, 00), the representation 

Dj (z; (s, 00)) = exp ( — / cr(x;z)dx) (2.23) 



in terms of the Jimbo-Miwa-Okamoto c-form of Painleve IV, namely 

°xx = 4(c r — xcr x ) 2 — 4o%(cr x + In), (2.24a) 

jn—\J2.n—7. 

cr(x;z) ~ z e~- T (x -> 00). (2.24b) 

As mentioned in the introduction, for the determinant (2.14b) and / = (0, s), Jimbo 
et al. (1980) found the representation (see also Tracy and Widom 1993a, Thm. 9) 

(bulk), , n - r gfoz) 



D™(z;(0,s))=exp^-y o dx j (2.25) 

in terms of the Jimbo-Miwa-Okamoto c-form of Painleve V, namely 

(xcr xx ) 2 = 4(cr — xcr x )(xa x — a — o^), (2.26a) 
z z 2 

c(x;z) ~ — iH ~x 2 (x — > 0). (2.26b) 

Finally there is Tracy and Widom's (1993b, 1994a) famous representation of the 
determinant (2.15b) for / = (s, 00), 

D^ ott \z; (s,oo)) = exp (— J (x — s)u(x;z) 2 dx^j (2.27) 

in terms of Painleve II, 5 namely 

u xx = 2« 3 + xu, u(x;z) ~ y/zAi(x) (x — >■ 00). (2.30) 

2.3.2. GOE and GSE in the bulk. With c(x) = cr(x;l) from (2.26) there holds the 
representation (Basor, Tracy and Widom 1992) 



E l( 0;s)=exp Ufy/^te Ws)"* (2.31a) 



E 4 (0;s/2)=cosh^y o y A ?M rf x j E 2 (0;s) 1/2 . (2.31b) 
Painleve representations for Ej (fc;s) and £4^; s) can be found in Basor et al. (1992). 



5 For a better comparison with the other examples we recall that Tracy and Widom (1994a, Eq. (1.16)) 
also gave the respresentation 



D2 Soft '(z; (s, 00)) = exp ^— J tr(x; z) dx J 



(2.28 



in terms of the Jimbo-Miwa-Okamoto a-foim of Painleve II: 

o* x = —4cr x {cr - xa x ) - 4o%, cr(x;z) ~ z (Ai'(x) 2 - xAi(x) 2 ) (x -> 00). (2.29) 
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2.3.3. GOE and GSE at the soft edge. Tracy and Widom (1996) found, with u(x) = 
w(x;l) being the Hastings-McLeod solution of (2.30), the representation 

Jl(l;s)=exp ^-^ ° U (x)dx^F 2 (l;s) 1/2 / (2.32a) 

F 4 (l;s) =coshQ^° U (x)dx^F 2 (l;s) 1/2 . (2.32b) 

More general, Dieng (2005) found Painleve representations of ¥\(k;s) and F^(k;s). 

2.4. Laguerre Ensembles. Here, we take, on x £ (0, 00) with parameter a > — 1, 
the weight functions 6 

• zv a (x) = x a e~ x ^ 2 for & = 1, the Laguerre Orthogonal Ensemble (LOE), 

• w a (x) = x a e~ x for & = 2, the Laguerre Unitary Ensemble (LUE), 

• w a {x) = x a e~ x for /3 = 4, the Laguerre Symplectic Ensemble (LSE). 
We define, for an open interval / C (0, 00), the basic quantity 

(n) 

iS-lev^' -^' a ) = ]P( exactr y k eigenvalues of the 

n x n Laguerre /3-ensemble with parameter a lie in /). (2.33) 

2.4.1. Scaling limits. The large matrix limit at the hard edge is (see Forrester 2006, 
p. 2993) 

f lim E < ?} v (k}4- 1 ir 1 J,x) B = 1 or B = 2, 

£ r»(w,«)= ~ 1 ; (.34) 

The large matrix limit at the soft edge gives exactly the same result (2.4) as for the 
Gaussian ensembles, namely (see Forrester 2006, p. 2992) 



^ = to Ej&(te4« + 2(2»)"*/,«) = 1 or „ = 2, ^ 

* ' |jimE^ ) (/c;4«+2(2n) 1/3 /, a ) B = 4, 

independently of a. Likewise, a proper bulk scaling limit yields Eg (A:; /) as for 
the Gaussian ensembles (see Tracy and Widom 1994b, p. 291). 

In the rest of this section we confine ourselves to the discussion of the LUE; 
determinantal formulae for LOE and LSE at the hard edge are given in Section 7. 

2.5. Determinantal Representations for the LUE. Here, the basic formula is (see 
Mehta 2004, §19.1) 



4ij E (^/^) = 1 #^D[ T ) E ( Z ;/, a ) 



(2.36a) 

z=l 



D LUe( Z ' /' a ) = det f 1 ~ zK n A f 1,2(7)) / (2.36b) 



6 We follow Forrester and Rains (2001) in this particular choice of the weights; as for the Gaussian 
ensembles notations differ from reference to reference by various scaling factors. 
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with the Laguerre kernel (the second form follows from Christoffel-Darboux) 

K nA (x,y) = ( x )n W = -Jn(n + oc) — — 

k=o x y 

(2-37) 

that is built from the Laguerre polynomials li*^ (x) by the L 2 (0, oo)-orthonormal 
systems of functions 



4 a) W = J r( fc + l + i) * a/2 ^ /2 4 a) (*)• (-38) 



The scaling limit at the hard edge is given by (Forrester 1993) 



(2.39a) 



2=1 



D< hard) (z; /, a) = det (l - zK a f L2(/) ) , (2.39b) 
with the Bessel kernel 

K a (x, y ) = ki^hMiVy) z^M^MVyl, (2 . 39c) 

Note that for non-integer parameter a the Laguerre kernel K„ A (x,y) and the 
Bessel kernel K a {x,y) exhibit algebraic singularities at x — or y — 0; see Sec- 
tion A.i for the bearing of this fact on the choice of numerical methods. 

2.6. Painleve Representations for the LUE. Tracy and Widom (1994c) calculated, 
for the determinant (2.36b) with / = (0, s), the representation 

Dg E (z; (0,s),«) = exp (- j° ^± dx^j (2.40) 

in terms of the Jimbo-Miwa-Okamoto c-form of Painleve V, namely 

{xctxx) 2 = {<r- xcr x - 2al + [In + a)cr x ) 2 - 4cr 2 (a x -n)(a x -n- a), (2.41a) 

r(n + a + l) K+1 , 
a{X ' Z) ~ Z r(n)r(« + l)r(« + 2) * {X 0) - (2 ' 4lb) 

Accordingly Tracy and Widom (1994b) obtained, for the determinant (2.39b) of 
the scaling limit at the hard edge, the representation 

Df ard) (z; (0,s), a ) = exp (- jf dx\ (2.42) 

in terms of the Jimbo-Miwa-Okamoto c-form of Painleve III, namely 

{xcr xx ) 2 = oc 2 cr 2 - cr x (cr - xa x )(Aa x - 1), (2.43a) 

a{x ' z) ~ T(a + 1)1 (a + 2) (I)^ ( * 0) ' { ^ h) 



ON THE NUMERICAL EVALUATION OF DISTRIBUTIONS IN RANDOM MATRIX THEORY 



li 



3. Numerics of Painleve Equations: the Need for Connection Formulae 

3.1. The Straightforward Approach: Solving the Initial Value Problem. All the 

five examples (2.24), (2.26), (2.30), (2.41), and (2.43) of a Painleve representation 
given in Section 2 take the form of an asymptotic initial value problem (IVP); that 
is, one looks, on a given interval (a,b), for the solution u(x) of a second order 
ordinary differential equation 

u"(x) = f(x,u(x),u'(x)) (3.1) 

subject to an asymptotic "initial" (i.e., one sided) condition, namely either 

u(x) ~ u a (x) (x — s- a) (3.2) 

or 

u{x)~u b (x) (x^b). (3.3) 
Although we have given only the first terms of an asymptotic expansion, further 
terms can be obtained by symbolic calculations. Hence, we can typically choose 
the order of approximation of u a {x) or U\,{x) at the given "initial" point. Now, the 
straightforward approach for a numerical solution would be to choose fl+ > a or 
b- < b sufficiently close and compute a solution v(x) of the initial value problem 

v"(x) = f{x,v(x),v'{x)) (3.4) 
subject to proper initial conditions, namely either 

v(a+) = u a (a+), v'{a+) = u' a (a+), (3.5) 

or 

v{b^) = u h (b_), v'(b-)=u' b (b^). (3.6) 
However, for principal reasons that we will discuss in this section, the straightfor- 
ward IVP approach unavoidably runs into instabilities. 

3.1.1. An example: the Tracy-Widom distribution F2. From a numerical point of 
view, the Painleve II problem (2.30) is certainly the most extensively studied case. 
We look at the Hastings-McLeod solution u(x) = u(x; 1) of Painleve II and the 
corresponding Tracy-Widom distribution 7 

F2(s)=exp^— J (x — s)u(x) 2 dx^j . (3.7) 

The initial value problem to be solved numerically is 

v"(x) =2v(xf + xv(x), v(b-) = Ai(fc-), »'(&_) = Ai'(fc-). (3.8) 

Any value of b- ^8 gives initial values that are good to machine precision (in 
IEEE double precision, which is about 16 significant decimal places). We have 
solved the initial value problem with b- — 12, using a Runge-Kutta method 
with automatic error and step size control as coded in Matlab's ode45, which is 
essentially the code published in Edelman and Persson (2005) and Edelman and 
Rao (2005). The red lines in Figure 1 show the absolute error \v(x) — u(x)\ and 
the corresponding error in the calculation of F?. We observe that the error of v(x) 



7 There is no need for a numerical quadrature here (and in likewise cases): simply add the differential 
equation (logics))" = — u(s) 2 to the system of differential equations to be solved numerically; see 
Edelman and Persson (2005). The same idea applies to the BVP approach in Section 3.2.1. 
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a. error in evaluating Fz(s) b. error in evaluating u(x) 

Figure i. Absolute error in evaluating the Tracy- Widom distribu- 
tion F2(s) and the Hastings-McLeod solution u(x) of Painleve II 
using different numerical methods; red: initial value solution (Mat- 
lab's ode45 as in Edelman and Persson 2005), which breaks down 
at about x = —5.56626; brown: boundary value solution (Matlab's 
bvp4c as in Dieng 2005), blue: boundary value solution by spec- 
tral collocation (Driscoll, Bornemann and Trefethen 2008); green: 
numerical evaluation of the Airy kernel Fredholm determinant 
(Bornemann 2010a), see also Section 4 (there is no u(x) here). The 
dashed line shows the tolerance 5 ■ 10 -15 used in the error control 
of the last two methods. All calculations were done in IEEE double 
precision hardware arithmetic. 



grows exponentially to the left of b- and the numerical solution ceases to exist 
(detecting a singularity) at about x = —5.56626. Though the implied values of F2 
are not completely inaccurate, there is a loss of more than 10 digits in absolute 
precision, which renders the straightforward approach numerically instable (that 
is, unreliable in fixed precision hardware arithmetic). 

To nevertheless obtain a solution that is accurate to 16 digits Prahofer and 
Spohn (2004) turned, instead of changing the method, to variable precision software 
arithmetic (using up to 1500 significant digits in Mathematica) and solved the initial 
value problem with b- = 200 and appropriately many terms of an asymptotic 
expansion Uj,(x). Prahofer (2003) put tables of u(x), Fz(s) and related quantities to 
the web, for arguments from —40 to 200 with a step size of 1/16. We have used 
these data as reference solutions in calculating the errors reported in Figure 1 and 
Table 1. 

3.1.2. Explaining the instability of the IVP. By reversibility of the differential equation, 
finite precision effects in evaluating the initial values at X = b- can be pulled back 
to a perturbation of the asymptotic condition u(x) ~ Ai(x) for 1^00. That is, 
even an exact integration of the ordinary differential equation would have to suffer 
from the result of this perturbation. Let us look at the specific perturbation 



u(x;6) ~ 6 ■ Ai(x) (x -> 00) with 6 = 1 + e. 



(3-9) 
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Table 1. Maximum absolute error and run time of the methods 
in Figure 1. The calculation was done for the 401 values of F2(s) 
from s = — 13 to s = 12 with step size 1/16. The I VP solution is 
only available for the 282 values from s = —5.5625 to s = 12. All 
calculations were done in hardware arithmetic. 



method reference max. error run time 



IVP/Matlab's ode45 


Edelman and Persson (2005) 


9.0 


10- 


-5 


11 sec 


BVP/Matlab's bvp4c 


Dieng (2005) 


1.5 


10- 


-10 


3.7 sec 


BVP/ spectral colloc. 


Driscoll et al. (2008) 


8.1 


10" 


-14 


1.3 sec 


Fredholm determinant 


Bornemann (2010a) 


2.0 


10" 


-15 


0.69 sec 




Figure 2. Sensitivity of Painleve II with the asymptotic condition 
u(x) ~ (1 + e)Ai(x) (x — > 00) for e m 0. The calculation was 
done with variable precision software arithmetic. Observe the 
dependence of the asymptotic behavior, for x — > —00, on the sign 
of e. 



The results are shown, for e = ±10 -8 and e = ±10~ 16 (which is already below the 
resolution of hardware arithmetic), in Figure 2 (see also Clarkson 2006, Figs. 11/ 12). 
Therefore, in hardware arithmetic, an error of order one in computing the Hastings- 
McLeod solution u(x) = u{x;\) from the IVP is unavoidable already somewhere 
before x m —12. 

This sensitive behavior can be fully explained by the connection formulae of 
Painleve II on the real axis, see Clarkson (2006, Thms. 9.1/2) and Fokas, Its, Kapaev 
and Novokshenov (2006, Thms. 10.1/ 2). Namely, for the Painleve II equation (2.30), 
the given asymptotic behavior u(x;8) ~ 9 ■ Ai(x), 9 > 0, as x — > 00 implies 
explicitly known asymptotic behavior "in the direction of x — » —00": 
(1) < 9 < 1: 

u(x;9) = d(0)|x|- 1/4 sin(!|x| 3/2 - |d(0) 2 log \x\ - (p{9)^j 

+ 0(|x|- 7/10 ) (x^-00), (3.10) 

with 

d(9) 2 = -7i" 1 log(l -9 2 ) (3.11a) 
<p{9) = |d(0) 2 log(2) +ar g r(l - {d{9) 2 ) - f; (3.11b) 
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Table 2. Blow-up points xq(1 + e) of u(x; 1 + e) with e > 0. 



e 




x (l +e) 


10" 


4 


-5.40049 30292 23929- • • 


10" 


8 


-8.011336780474602- • • 


10" 


•12 


-10.21585052253541 • • • 


10- 


-16 


-12.1916 5664375788- •• 



(2) = 1: _ 

u(x;l) = ^J-^+0(x- 5/1 ) (*->-«>); (3.12) 

(3) 8 > 1: there is a pole at a finite Xq(9) 6 R, such that 

1 

— WflY (*4*o(0))- (3-i3) 

We observe that, for x — > — 00, the asymptotic behavior of the Hastings-McLeod 
solution w(x; 1) separates two completely different regimes: an oscillatory (0 < 1) 
from a blow-up solution (0 > 1). The blow-up points Xq(8) are close to the range 
of values of x which are of interest in the application to RMT, see Table 2. 

3.1.3. Explaining the separation of asymptotic regimes. The deeper reason for this 
separation property comes from comparing (2.27) with (2.15b), that is, from the 
equality 

det(J - 2 K Ai \ L 2( S)Bo) ) = exp f-J^ (x - s)u(x;8) 2 dx^j . (3.14) 

Here, Il 2 (s,oo) is a positive self-adjoint trace class operator with spectral radius 
(see Tracy and Widom 1994a) 

p(s) = A max (KAit L 2( S/ c°)) < 1 - ( 3-i5) 
Obviously, there holds p(s) — >■ for s — >■ 00. On the other hand, since for 9 = 1 the 
determinant (3.14) becomes the Tracy- Widom distribution F2(s) with F2( s ) ~^ as 
s — > —00, we conclude that p(s) — >■ 1 as s — >■ —00. 

Now, we observe that the determinant (3.14) becomes zero if and only if u(x; 6) 
blows up at the point x — s. By the Painleve property, such a singularity must be a 
pole. By Lidskii's theorem (Simon 2005, Thm. 3.7), the determinant (3.14) becomes 
zero if and only if 6~ 2 is an eigenvalue of K Ai \ L 2^ S/O0 y Therefore, a blow-up point 
of u(x;6) at x — s implies necessarily that 

e ^ P {s)- 1/2 > 1. (3.16) 

On the other hand, if 8 > 1 there must be, by continuity, a largest point s = x${8) 
such that 8 = p{s)~ 1 ^ 2 , which gives us the position of the pole of the connection 
formula. This way, using the methodology of Section 4, we have computed the 
numbers shown in Table 2. 

A similar line of arguments shows that in the other cases (2.24), (2.26), (2.41), 
and (2.43) of a Painleve representation given in Section 2, the parameter z = 1 
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(which is the most significant choice for an application to RMT) also belongs to a 
connecting orbit a(x; 1) that separates different asymptotic regimes. In particular, 
we get poles at finite positions if and only if z > 1. Hence, the numerical difficulties 
observed with the initial value approach have to be expected in general. 

3.2. The Stable Approach: Solving the Boundary Value Problem. The stable nu- 
merical solution of a connecting orbit separating different asymptotic regimes has 
to be addressed as a two-point boundary value problem (BVP), see, e.g., Deuflhard 
and Bornemann (2002, Chap. 8). That is, we use the information from a connection 
formula to infer the asymptotic for x — » b from that of x — »■ a, or vice versa, and to 
approximate u(x) by solving the BVP 

v"(x) = f(x,v(x),v'(x)), v(a+) = ««(«+), v(b-) = u h (b-). (3.17) 

Thus, four particular choices have to be made: The values of the finite boundary 
points fl+ and b-, and the truncation indices of the asymptotic expansions at x — »■ a 
and x — >■ b that give the boundary functions u a (x) and u\,{x). All this has to be 
balanced for the accuracy and efficiency of the final method. 8 

3.2.1. An Example: the Tracy-Widom distribution ^{x). Let us look, once more, at 
the Hastings-McLeod solution u(x) = u(x;l) of (2.30) and the corresponding 
Tracy- Widom distribution (3.7). By definition, we have 

u(x)~Ai(x) (x — > 00). (3-18) 

The asymptotic result for x — > —00 as given in the connection formula (3.12) is not 
accurate enough to allow a sufficiently large point a+ to be used. However, using 
symbolic calculations it is straightforward to obtain, from this seed, the asymptotic 
expansion (Tracy and Widom 1994a, Eq. (4.1)) 

ric / 1 _ 3 73 _ 6 10657 _ 9 13912277 _ 12 n . _ 15 ,\ 

U(x) = j- - 1 + -X - — X b +^-—rX 9 — - X 12 + 0(X 15 ) 

w V 2 V 8 128 1024 32768 v ' ) 

(x^-00). (3.19) 

Dieng (2005) chooses these terms as u a (x), as well as a + = — 10, = 6 and 
Mj,(x) = Ai(x). Using Matlab's fixed-order collocation method bvp4c, he calculates 
solutions whose errors are assessed in Figure 1 and Table 1. The accuracy is still 
somewhat limited and he reports (p. 88) on difficulties in obtaining a starting iterate 
for the underlying nonlinear solver. A more promising and efficient approach to 
obtain near machine precision is the use of spectral collocation methods. Because 
of analyticity, the convergence will be exponentially fast. This can be most elegantly 
expressed in the newly developed chebop system of Driscoll et al. (2008), a Matlab 
extension for the automatic solution of differential equations by spectral collocation. 
In fact, the evaluation of the Tracy- Widom distribution is Example 6.2 in that paper. 
Here, the first four terms of (3.19) are chosen as u a (x), as well as a + = —30, b- = 8 

Q 

An early variant of this connection-formula based approach can be traced back to the work of 
Wu, McCoy, Tracy and Barouch (1976, App. A): there, a Painleve III representation of the spin-spin 
correlation function of the two-dimensional Ising model was evaluated by joining a forward integration 
of the IVP from a + to some interior point c £ (a + , b- ) with a backward integration of the IVP from 
fc_ to c. The difference of the two IVP solutions at c was used as a rough error estimate. Though not 
quite a BVP solution, it is close in spirit. Actually, this was the approach originally used by Tracy and 
Widom (1993b, 1994a) in their numerical evaluation of F2 (personal communication by Craig Tracy). 
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and U),(x) = Ai(x). The Newton iteration is started from a simple affine function 
satisfying the boundary conditions; see Figure 1 and Table 1 for a comparison of 
the accuracy and run time. 

3.3. A List of Connection Formulae. For the sake of completeness we collect the 
connection formulae for the other Painleve representations (2.24), (2.26), (2.41), and 
(2.43). References are given to the place where we have found each formula; we 
did not try to locate the historically first source, whatsoever. Note that a rigorous 
derivation of a connection formula relies on deep and involved analytic arguments 
and calculations; a systematic approach is based on Riemann-Hilbert problems, 
see Deift (1999b) and Fokas et al. (2006) for worked out examples. 

• The Painleve III representation (2.43), for LUE with parameter oc at the 
hard edge, satisfies (Tracy and Widom 1994b, Eq. (3.1)) 

tr(*;l) = --|Vx + 0(l) (X-+00). (3-20) 

• The Painleve IV representation (2.24), for ^-dimensional GUE, satisfies 
(Tracy and Widom 1994c, Eq. (5.17)) 

o~(x;l) — — 2nx — nx' 1 + 0(x~ 3 ) (x — > — 00). (3-21) 

It is mentioned there that c(x;z) has, for z > 1, poles at finite positions. 
This is consistent with the line of arguments that we gave in Section 3.1.3. 

• The Painleve V representation (2.26), for the bulk scaling limit of GUE, 
satisfies (Basor et al. 1992, p. 6) 

/ 2 

x 1 

-r Z = 1 

a(x;z) ~ I 4 [x _> oo ) ( 3 . 22 ) 

-log(l-z)- 0<z<l 

• The Painleve V representation (2.41), for n-dimensional LUE with parame- 
ter a, satisfies (Forrester and Witte 2002, Eq. (1.42)) 

c(x; 1) = nx — an + an 2 x~ l + 0(x~ 2 ) (x — > 00). (3-23) 

3.4. Summary. Let us summarize the steps that are necessary for the numerical 
evaluation of a distribution function from RMT given by a Painleve representation 
on the interval (a, b) (that is, the second order differential equation is given together 
with an asymptotic expansion of its solution at just one of the endpoints a or b): 

(1) Derive (or locate) the corresponding connection formula that gives the 
asymptotic expansion at the other end point. This requires considerable 
analytic skills or, at least, a broad knowledge of the literature. 

(2) Choose a + > a and b_ < b together with indices of truncation of the 
asymptotic expansions such that the expansions themselves are sufficiently 
accurate in (a, a + ) and (£>_, b) and the two-point boundary value problem 
(3.17) can be solved efficiently. This balancing of parameters requires a 
considerable amount of experimentation to be successful. 

(3) The issues of solving the boundary value problem (3.17) have to be ad- 
dressed: starting values for the Newton iteration, the discretization of 
the differential equation, automatic step size control etc. This requires a 
considerable amount of experience in numerical analysis. 
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Thus, much work has still to be done to make all this a "black-box" approach. 

4. Numerics of Fredholm Determinants and Their Derivatives 

4.1. The Basic Method. Bornemann (2010a) has recently shown that there is an 
extremely simple, accurate, and general direct numerical method for evaluating 
Fredholm determinants. By taking an wz-point quadrature rule 9 of order 10 m with 
nodes Xj £ (a, b) and positive weights Ws, written in the form 

m 

E w jf( x j) s 

the Fredholm determinant 

d(z) = det(J 

is simply approximated by the corresponding wz-dimensional determinant 

d m (z) = det (Sij-zw} /z K(xi,Xj)wy 2 ^ . (4.3) 

This algorithm can straightforwardly be implemented in a few lines. It just needs 
to call the kernel K(x,y) for evaluation and has only one method parameter, the 
approximation dimension m. 

If the kernel function K(x,y) is analytic in a complex neighborhood of (a,b), 
one can prove exponential convergence (Bornemann 2010a, Thm. 6.2): there is a 
constant p > 1 (depending on the domain of analyticity of K) such that 

d m (z) - d(z) = 0(p- m ) (m 00), (4.4) 

locally uniform in 2 £ C. (Note that d(z) is an entire function and d m (z) a polyno- 
mial.) This means, in practice, that doubling m will double the number of correct 
digits; machine precision of about 16 digits is then typically obtained for a rel- 
atively small dimension m between 10 and 100. This way the evaluation of the 
Tracy-Widom distribution f2( s )/ at a given argument s, takes just a few millisec- 
onds; see Figure 1 and Table 1 for a comparison of the accuracy and run time with 
the evaluation of the Painleve representation. 

4.2. Numerical Evaluation of Finite-Dimensional Determinants. Let us write 

d m {z) = det(I - zA m ), A m £ K mxm , (4.5) 

for the finite-dimensional determinant (4.3). Depending on whether we need its 
value for just one z (typically z = 1 in the context of RMT) or for several values of 
z (such as for the calculation of derivatives), we actually proceed as follows: 

(1) The value d m {z) at a given point z £ C is calculated from the LU decompo- 
sition of the matrix I — zA m (with partial pivoting). Modulo the proper sign 
(obtained from the pivoting sequence), the value is given by the product of 
the diagonal entries of U (Stewart 1998, p. 176). The computational cost is 
of order 0(m 3 ), including the cost for obtaining the weights and nodes of 
the quadrature method, see Bornemann (2010a, Footnote 5). 



f(x) dx, 

zK \l 2 (a,b)) 



(4-i) 



(4-2) 



9We choose Clenshaw-Curtis quadrature, with a suitable meromorphic transformation for (semi) 
infinite intervals (Bornemann 2010a, Eq. (7.5)). For the use of Gauss-Jacobi quadrature see Section A.i. 
I0 A quadrature rule is of order m if it is exact for polynomials of degree m — 1. 
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(2) The polynomial function d m (z) itself is represented by 

m 

d m (z)=Yl(l-z\ j (A m )). (4.6) 

;=i 

Here, we first calculate the eigenvalues Aj(A m ) (which is slightly more 
expensive than the LU decomposition, although the computational cost of, 
e.g., the QR algorithm is of order 0(m 3 ), too). The subsequent evaluation 
of d m (z) costs just 0(m) operations for each point z that we care to address. 

4.3. Numerical Evaluation of Higher Derivatives. The numerical evaluation of 
expressions such as (2.11a) requires the computation of derivatives of the determi- 
nant d(z) with respect to z. We observe that, by well known results from complex 
analysis, these derivatives enjoy the same kind of convergence as in (4.4), 

d%\z)-dW(z)=0(p- m ) (m ->«>), (4.7) 

locally uniform in z £ C, with an arbitrary but fixed k = 0, 1, 2, . . . 

The numerical evaluation of higher derivatives is, in general, a badly conditioned 
problem. However, for entire functions / such as our determinants we can make 
use of the Cauchy integrals 11 

f {k) (z) = ^ e- ike f(z + re ie ) d9 (r > 0). (4.8) 

Since the integrand is analytic and periodic, the simple trapezoidal rule is expo- 
nentially convergent (Davis and Rabinowitz 1984, §4.6.5); that is, once again, p 
quadrature points give an error O(p^P) for some constant p > 1. 

Theoretically, all radii r > are equivalent. Numerically, one has to be very 
careful in choosing a proper radius r (for a detailed study see Bornemann 2009). 
The quantity of interest in controlling this choice is the condition number of the 
integral, that is, the ratio 



f n e - ike f(z + re w ) d6 / [ " \f(z + re w )\ d6 . (4.9) 
Jo Jo 

For reasons of numerical stability, we should choose r such that k ps 1. Some 
experimentation has led us to the choices r = 1 for the bulk and r — 0.1 for the 
edge scaling limits (but see also Example 12.3 in Bornemann 2009). 

4.4. Error Control. Exponentially convergent sequences allow us to control the 
error of approximation in a very simple fashion. Let us consider a sequence d m — » d 
with the convergence estimate 

d m -d = o( P - m ) (4.10) 

for some constant p > 1. If the estimate is sharp, it implies the quadratic conver- 
gence of the contracted sequence d^, namely 

|d 2 , ;+ i — d\ < c|d 2 <; — d\ 2 (4-n) 



lx For more general analytic / one would have to bound the size of the radius r to not leave the 
domain of analyticity. In particular, when evaluating (2.16) we have to take care of the condition 
r < nun | (1 — A) /A|, where A runs through the eigenvalues of the matrix kernel operator. 
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Table 3. Approximation of the Airy kernel determinant F2(— 2) = 
det(7 — ^Aiti 2 (-2,oo)) by (4.3), using m-point Clenshaw-Curtis 
quadrature meromorphically transformed to the interval (—2, 00), 
see Bornemann (2010a, Eq. (7.5)). Observe the apparent quadratic 
convergence: the number of correct digits doubles from step to 
step. Thus, in exact arithmetic, the value for m = 64 would be 
correct to about 20 digits; here, the error saturates at the level of 
machine precision (2.22 • 10~ 16 ): all 15 digits shown for the m = 64 
approximation are correct. 



m 


dm 


\d m — £2 


("2 


)l 


error estimate (4.12) 


8 


0.38643 7295515158 


2.67868 ■ 


10- 


2 


2.67817 ■ 10- 2 


16 


0.413219001146910 


5.14136 • 


10- 


6 


5.14138 • 10- 6 


32 


0.413224142527728 


2.26050 ■ 


10- 


11 


2.26046 ■ 1Q- 11 


64 


0.413224142505123 


4.44089 ■ 


10- 


16 





for some c > 0. A simple application of the triangle inequality gives then 

\d 2 i - d\ < 1 / 2q+1 \, - - rf 2'/+! I (l °°)- (4- 12 ) 
1 — c\Q-m — a 1 

Thus we take | d2i — d 2 q+i | as an excellent error estimate of \d2i — d\ and as a quite 
"conservative" but absolutely reliable estimate of \d 2 q+i — d\. Table 3 exemplifies 
this strategy for the calculation of the value F2{— 2). 

4.5. Numerical Evaluation of Densities. The numerical evaluation of the proba- 
bility densities belonging to the cumulative distribution functions F(s) given by a 
determinantal expression requires a low order differentiation with respect to the 
real-valued variable s (which cannot easily be extended numerically into the complex 
domain). Nevertheless these functions are typically real-analytic and therefore 
amenable to an excellent approximation by interpolation in Chebyshev points. To 
be specific, if F(s) is given on the finite interval [a, b] and So, . . . ,s m denote the 
Chebyshev points of that interval, the polynomial interpolant p m (s) of degree m is 
given by Salzer's (1972) barycentric formula 

E^ (-l) /C f(Sic)/(s- S/0 

where the double primes denote trapezoidal sums, i.e., the first and last term 
of the sums get a weight 1/2. This formula enjoys perfect numerical stability 
(Higham 2004). If F is real analytic, we have exponential convergence once more, 
that is 

||F- Pm||oo =0{p- m ) (m^oo) (4.14) 

for some constant p > 1 (see Berrut and Trefethen 2004). Low order derivatives 
(such as densities) and integrals (such as moments) can easily be calculated 
from this interpolant. All that is most conveniently implemented in Battles and 
Trefethen's (2004) chebf un package for Matlab (see also Driscoll et al. 2008). 



«oo = ^t ,; m ,\:, : > (4-13) 
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a. smallest eigenvalue of LUE (n = 80, a — 40) b. largest eigenvalue of LUE (n = 80, a = 40) 



Figure 3 . CDF (green) and PDF (dark blue) of the smallest and 
largest eigenvalue for n-dimensional LUE with parameter a with 
n = 80, a = 40. Also shown are the scaling limits at the hard and 
soft edge (CDF in red, PDF in light blue). 

4.6. Examples. We illustrate the method with three examples. More about the 
software that we have written can be found in Section 9. 

4.6.1. Distribution of smallest and largest level in a specific LUE. We evaluate the 
cumulative distribution functions (CDF) and probability density functions (PDF) 
of the smallest and largest eigenvalue of the n-dimensional LUE with parameter a 
for the specific choices n = 80 and a. = 40. (Note that for parameters of this size 
the numerical evaluation of the Painleve representation (2.41) becomes extremely 
challenging.) Specifically, we evaluate the CDFs (the PDFs are their derivatives) 

P(A min < s) = 1 - Effi E (0; (0,s),«) (4.15) 

and 

P(A ma x < 4n + 2a + 2 + 2(2n) 1/3 s) 

= Eg E (0;(4n + 2a+2 + 2(2n) 1/3 s,oo),a). (4.16) 

Additionally, we calculate the CDFs of the scaling limits; that is, 

1-E^ hard) (0;(0,4ns), a ) (4.17) 

at the hard edge and the Tracy-Widom distribution 

F 2 (s)=Ef ft) (0;(s,co)) (4.18) 

at the soft edge. All that has to be done to apply our method is to simply code the 
Laguerre, Bessel, and Airy kernels. Figure 3 visualizes the functions and Table 4 
shows their moments to 10 correct decimal places. 
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Table 4. Moments of the distributions (4.15), (4.16), (4.17), and 
(4.18) for LUE with n = 80 and a = 40. We show ten correctly trun- 
cated digits (that passed the error control) and give the computing 
time. All calculations were done in hardware arithmetic. 



CFD 


mean 


variance 


skewness 


kurtosis 


time 


(4-15) 


5.14156 81318 


0.3434752478 


0.0431330951 


-0.0292563564 


1.0 sec 


(4-17) 


6.3558698372 


0.5210615307 


0.04102 67718 


-0.0294322640 


1.2 sec 


(4.16) 


-2.4391384563 


0.8934123428 


0.26271 64962 


0.1278351672 


4.1 sec 


(4.18) 


-1.7710868074 


0.8131947928 


0.2240842036 


0.0934480876 


1.0 sec 




a. ir-level spacing functions E 2 (k; s) of GUE b. density d s F 2 {k;s) of fc-th largest level in GUE 

Figure 4. Plots of the A:-level spacing functions E2 {k; s) in the bulk 
scaling limit of GUE (k = 0, . . . , 14; larger k go to the right), and 
of the probability density functions d s F2(k;s) of the A:-th largest 
level in the soft edge scaling limit of GUE (k — 1, ... ,6; larger 
k go to the left). The underlying calculations were all done in 
hardware arithmetic and are accurate to an absolute error of about 
5 ■ 10~ 15 . Each of the two plots took a run time of about 30 seconds. 
Compare with Mehta (2004, Fig. 6.4) and Tracy and Widom (1994a, 
Fig. 2). 



4.6.2. The distribution of k-level spacings in the bulk of GUE. By (2.5) and (2.14), 
the fc-level spacing functions in the bulk scaling limit of GUE are given by the 
determinantal expressions 



E 2 (k;s) 



(-1)* 



k\ dz k 



det \l-zK. 



sinlL 2 (0,s) 



(4-!9) 



z=l 



Mehta and des Cloizeaux (1972) evaluated them using Gaudin's method (which 
be briefly described in Section 1.1.1); a plot of these functions, for k from o up 
to 14, can also be found in Mehta (2004, Fig. 6.4). Now, the numerical evaluation 
of the expression (4.19) is directly amenable to our approach. The results of our 
calculations are shown in Figure 4.a. 
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4.6.3. The distributions of the k-th largest level at the soft edge ofGUE. By (2.9) and 
(2.15), the cumulative distribution functions of the A:-th largest level in the soft 
edge scaling limit of GUE are given by the determinantal expressions 



fc-i f—iv di / 

F 2( fc ; g ) = E -jr^j det i. 1 -^Aib (s ,oo) ) ^, (4.20) 



which is directly amenable to be evaluated by our numerical approach. The results 
of our calculations are shown in Figure 4-b. 

Remark. To our knowledge, prior to this work, only calculations of the particu- 
lar cases k = 1 (the largest level) and k — 2 (the next-largest level) have been 
reported in the literature (Tracy and Widom 2000, Dieng 2005). These calcula- 
tions were based on the representation (2.27) of the determinant in terms of the 
Painleve II equation (2.30). The evaluation of F2( s ) = ft(l;s) was obtained from 
the Hastings-McLeod solution u(x) = u(x; 1), see Section 3.2.1. On the other hand, 
the evaluation of F2(2;s) required the function 

zv(x) = d z u(x;z)\ z= i, (4.21) 

which, by differentiating (2.30), is easily seen to satisfy the linear ordinary differ- 
ential equation 

w"(x) = (6u(x) 2 + x)w(x), w{x) ~ jAi(x) (x — > 00). (4.22) 

Obtaining the analogue of the connection formula (3.19) requires some work 
(though, since the underlying differential equation is linear, it poses no fundamen- 
tal difficulty) and one gets (see Tracy and Widom (1994a, p. 164) who also give 
expansions for larger k) 

e2 ^ ( ~* )3/2/3 (, I 7 / \-3/2 1513 _ 3 ^-9/2^ 

(x->-co). (4.23) 

Note that the exponential growth points, once more, to the instability we have 
discussed in Section 3.1.2. 



5. The Distribution of k-Level Spacings in the Bulk: GOE and GSE 

Mehta (2004, Chap. 20) gives determinantal formulae for the A:-level spacing 
functions Ep(k;s) in the bulk scaling limit that are (also in the cases /S = 1 and 
f> = 4 of the GOE and GSE, respectively) directly amenable to the numerical 
approach of Section 4. These formulae are based on a factorization of the sine 
kernel determinant (2.14b), which we describe first. 

Since K s { n is a convolution operator we have the shift invariance 

det (l - zJC s iJ L 2 (0/2 f)) = det (l - zK sin \ L2{ _ t /f) ) . (5.1) 

Next, there is the orthogonal decomposition L 2 (-t,t) = Xf ven © Xf dd into the 
even and odd functions. On the level of operators, this corresponds to the block 
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diagonalization 

K sm\ L 2(-t,t)= ( Vxf ve »©X° dd (5- 2 ) 

\ sin/ 

with the kernels 

K in( x 'V) = l{Ksm{x,y)±K sin (x,-y)). (5.3) 
Further, there is obviously 

x+r -( KL V irr -f° 

^sin* L 2 (-t,t)~ I I I Xf ven ©X° dd ' ^sin I L 2 (-t,t)~ I ^_ I I Xf ven ffiX° dd • 

(54) 

Hence, we get the factorization 

det(l-zK sln r L2( _ M) ) =det(l-zK+ n r L2( _ M) ) det (i - zK^ \ L 2 { _ tit) ) . (5.5) 
Now, upon introducing the functions 



E ± (k;s) = Ul!^det(l-z^J L2( _ s/2 , /2) ) 



(5-6) 



z=l 



we obtain, using the factorization (5.5) and the Leibniz formula applied to (4.19), 
the representation 

k 

E 2 (k;s) = £E + (/;s)E_(*-;;s) (5.7) 

j=0 

of the A:-level spacing functions in the bulk of GUE. The important point here is 
that Mehta (2004, Eqs. (20.1.20/21)) succeeded in representing the A:-level spacing 
functions of GOE and GSE in terms of the functions E±, too: 

E 1 (0;s) = E+(0;s), (5.8a) 

Ei(2fc-l;s) = E_(Jfc-l;s)-E 1 (2Jfc-2;s), (5.8b) 

E 1 (2fc;s) = E + {k;s) - E r (2k - l;s), (5.8c) 

(k = 1,2,3,...) for GOE and 

E 4 (k;s) = i(£ + (fc;2s)+E_(fc;2s)) (5.9) 

(k — 0, 2, 3, . . .) for GSE. Based on these formulae, we have used the numerical 
methods of Section 4 to reproduce the plots of Mehta (2004, Figs. 7.3/11.1). The 
results of our calculations are shown in Figure 5. 

6. The k-th Largest Level at the Soft Edge: GOE and GSE 

In this section we derive new determinantal formulae for the cumulative dis- 
tribution functions Fj(fc;s) and F±(k;s) of the fc-th largest level in the soft edge 
scaling limit of GOE and GSE. We recall from (2.9) that 

k-\ 

Fp(k;s) = s), (6.1) 

j=0 

where we briefly write 

E^;s)=E<j Soft) (fc;( S ,oo)). (6.2) 
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2 4 6 8 10 12 14 2 4 6 8 10 12 14 

s S 



a. fc-level spacing functions Ej(/c;s) of GOE b. /:-level spacing functions Ei(k;s) of GSE 

Figure 5. Plots of the A:-level spacing functions Ei(fc;s) in the 
bulk scaling limit of GOE and E^fas) in the bulk scaling limit 
of GSE (k = 0, ... ,14; larger k go to the right). The underlying 
calculations were all done in hardware arithmetic and are accurate 
to an absolute error of about 5 • 10~ 15 . Each of the two plots took 
a run time of less than one minute. Compare with Mehta (2004, 
Figs. 7.3/11.1). 



The new determinantal formulae of this section are amenable to the efficient 
numerical evaluation by the methods of Section 4; but, more important, they 
are derived from a determinantal equation (6.10) whose truth we established by 
numerical experiments before proving it rigorously. Therefore, we understand this 
section as an invitation to the area of Experimental Mathematics (Borwein and 
Bailey 2004). 

In a broad analogy to the previous section we start with a factorization of the 
Airy kernel determinant (2.15b) that we have learnt from Ferrari and Spohn (2005, 
Eq. (34)). Namely, the integral representation 

/■CO 

K Ai (x,y) = Ai(jc + £)Ai(y + £)d£ (6.3) 
Jo 

of the Airy kernel implies, by introducing the kernels 

T s {x,y) = Ai(x + y + s), V M (x,y) = ^Ai ' (6 ' 4) 

the factorization 




= det (j - T s r L2(0/Oo) ) ■ det (i + T s t L 2 (aoo) ) 

= det (i - v^V Ai r L2(s co) ) ■ det (i + V A ib( s ,oo)) (6-5) 
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that is valid for the complex cut plane z 6 C \ (— oo,0]. Now, upon introducing the 
functions 

£±(k;s) = { ~^^ k det(lT^~zV M \ LHSi0o) ) (6.6) 

z=l 

we obtain, using the factorization (6.5) and the Leibniz formula applied to (2.15), 
the representation 

k 

£ 2 (*;s) = ££+(/;*)£_(* -;;s). (6.7) 
;=0 

Further, Ferrari and Spohn (2005, Eqs. (33/35)) proved that 

£i(0;s) = E+(0;s) (6.8a) 

£ 4 (0;s) = i(£+(0;s) + £_(0;s)). (6.8b) 

The similarity of the pairs of formulae (6.7)7(5.7), (6.8a)/ (5.8a), and (6.8^/(5.9) 
(the last with k = 0) is absolutely striking. So we asked ourselves whether (6.8b) 
generalizes to the analogue of (5.9) for general k, that is, whether 

E 4 (*;s) = s)+E_(*;s)) (jfc = 0,1,2,. . .) (6.9) 

is valid in general. In view of (2.16) such a result is equivalent to the following 
theorem. We first convinced ourselves of its truth in the sense of experimental 
mathematics: by numerically 12 checking its assertion for 100000 randomly chosen 
arguments. Thus being encouraged, we then worked out the proof given below. 13 

Theorem 6.1. The determinantal equation 
D 4 (z; (s,~)) 1/2 

= \ (det(l-^V Ai f L2(S/Co) ) +det(l + ^V Ai r i 2 (sH )) (6.10) 

holds for all s 6 R and z in the complex domain of analyticity that contains z — 1. 

Proof. Since the operator theoretic arguments of Ferrari and Spohn (2005) cannot 
directly be extended to yield (6.10) we proceed by using Painleve representations. 
Dieng (2005, Thm. 1.2.1, Eq. (1.2.2)) proved that (2.32b) generalizes to 

D 4 (0 2 ;(s,co)) 1 / 2 = D l f oh) {9 2 ;{s,oo)) l ' 2 cosh(^- u{x;9)dx^j (6.11a) 

in terms of the Painleve representation 

u xx =2u 3 +xu, u{x;6) ~ 6 ■ Ai(x) (x — > 00). (6.11b) 

Here we put y/z — 9 and observe that (6.11a) obviously extends, by the symmetry 
of the Painleve II equation, from < 9 1 to the range — 1 $C 9 1. In view 



I2 The function 04(2; (s, 00)) was evaluated using the extension of our method to matrix kernel 
determinants that will be discussed in Section 8.1; see also Example 9.2.2 for a concrete instance. 
13 Later though, we found that the result has recently been established by Forrester (2006, Eq. (1.23)). 
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of (2.28) Dieng's representation (6.11a) readily implies, by analytic continuation, 
the asserted formula (6.10) if the representation 

det (l - 8 VAibf^co)) = exp (~\ £ u(x;8)dx^j det (j - 8 2 K Ai \ L 2 {Si0o) ) 1/2 

= exp ^ — (m(x; #) + (x — s)m(x; #) 2 ) dx^j (6.12) 

happens to be true for all —1 1 (note that we have chosen the signs in 
accordance with the the special cases 8 = ±1 as given by (6.8) and (2.32)). Now, 
this particular Painleve representation can directly be read off from the work of 
Desrosiers and Forrester (2006, Eqs. (4.8/19)), which completes the proof. □ 

It remains to establish formulae for the GOE functions E\(k;s) that are struc- 
turally similar to (5.8). To this end we use the interrelationships between GOE, GUE, 
and GSE found by Forrester and Rains (2001, Thm. 5. 2), 14 which can symbolically 
be written in the form 

GSE„ = even(GOE 2 „+i), (6.13a) 

GUE„ = even(GOE„ U GOE„ +1 ). (6.13b) 

The meaning is as follows: First, the statistics of the ordered eigenvalues of the 
n-dimensional GSE is the same as that of the even numbered ordered eigenvalues 
of the In + 1 -dimensional GOE. Second, the statistics of the ordered eigenvalues 
of the n-dimensional GUE is the same as that of the even numbered ordered 
levels obtained from joining the eigenvalues of a n-dimensional GOE with the 
eigenvalues of a statistically independent n + 1 -dimensional GOE. 

Now, (6.13a) readily implies, in the soft edge scaling limit (2.4), that the cumu- 
lative distribution function of the A:-th largest eigenvalue in GSE agrees with the 
cumulative distribution function of the 2A:-th largest of GOE, 

F i (k;s)=F l (2k;s) (k = 1,2,3, .. .), (6.14) 

the so-called interlacing property. Therefore, in view of (6.1) and (6.9) we get 

E l (2k;s)+E 1 (2k + l;s) = \{E + {k;s) + £_(*;s)). (6.15) 

Further, the combinatorics of (6.13b) implies, in the soft edge scaling limit (2.4): 
exactly k levels of GUE are larger than s if and only if exactly 2k or 2k + 1 levels 
of the union of GOE with itself are larger than s. Here, j levels are from the first 
copy of GOE and 2k — j, or 2k + 1 — j, are from the second copy. Since all of these 
events are mutually exclusive, we get 

2k 2k+l 
E 2 (k;s) = '£ j E 1 (j;s)E 1 (2k-j;s)+ £ E 1 (;';s)£ 1 (2fc + 1 - j;s). (6.16) 
;'=0 j=0 

Finally, the following theorem gives the desired (recursive) formulae for the 
functions Ej(fc;s) in terms of E±(k;s). Note that these recursion formulae, though 
being quite different from (5.8), share the separation into even and odd numbered 
cases. 



I 4That is why we have chosen, in defining the Gaussian ensembles, the same variances of the 
Gaussian weights as Forrester and Rains (2001). 
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Theorem 6.2. The system (6.y), (6.8a), (6.15), and (6.16) of functional equations has the 
unique, recursively defined solution 

( 2 !\ 



• ' q) 

7=0 



E a (2k; s) = E + (k; s) - £ g^T^y Ei (2* - 2; ~ 1; s), (6.17a) 



£1 (2* + 1; s) = £+(fc;S) + £ - (fc;S) - Ei (2k; s). (6.17b) 

Prao/ We introduce the generating functions 

00 

feven(x) = £ Ei (2fc; s)x 2,c (6.18a) 

00 

/odd(^) = E £ i ( 2ic + 1; ( 6 - l8b ) 

fc=0 
co 

?± W = [E±(t; S )r a (6.18c) 

Ferrari and Spohn's (2005) representation (6.8a) translates into the constant term 
equality (which breaks the symmetry of the other functional equations) 

/even(O) =g+(0). (6.1 9 ) 

Equating (6.7) and (6.16) translates into 

/even(^) 2 +2/even(x)/ odd (x)+X 2 /odd(^) 2 =£+(*)•£-(*)• (6-20) 

Finally (6.15) translates into 

/even (*) + / od d (*) = \ (g+ (x) + g- (x)). (6.2l) 
Elimination of g- from the last two equations results in the quadratic equation 

(feven(x) - g+(x)) 2 +2(f even (x) - g + (x))f odd (x) + X 2 f odd {x) 2 = 0. (6.22) 

Solving for f eV en(x) — g+ (x) gives the two possible solutions 



feven(x) ~ g+(x) = - ( 1 ± \J \ - X 2 ) f odd (x). (6.23) 



Because of / o dd(0) = Ei(l;s) > we have to choose the negative sign of the 
square root to satisfy (6.19). To summarize, we have obtained the mutual relations 

feven(x) = g+ (x) - (l - \J \ - X 2 ) f odd (x), (6.24a) 
/oddM = \ (g+(x) +g-(x)) -feven(x), (6.24b) 

which then, by the Maclaurin expansion 

< 2 h 

1 - vT 3 *? = £ ^FiyTT)^'"' < 6 ^ 

translate back into the asserted recursion formulae of the theorem. □ 
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a. density d s F\(k;s) of k-th largest level in GOE b. absolute error of Dieng's (2005) evaluation 



Figure 6. Left: plots of the probability densities of the A:-th largest 
level in the soft edge scaling limit of GOE (k = 1, . . . , 6; larger k go 
to the left) based on the recursion formulae (6.17). The underlying 
calculations were all done in hardware arithmetic and are accurate 
to an absolute error of about 5 • 10~ 15 (dashed line). Right: the 
absolute error (taking the values of the calculations on the left as 
reference) of the Painleve II based calculations by Dieng (2005). 



Remark 6.3. The system (6.24) can readily be solved for f eV en{x) and f dd( x ) to 
yield 

00 

E 1 (k;s)x k = feven(x) + xf odd (x) 

k=0 

= I (*+M f 1 + VttI) + z- (x) ( ] - \/ttI)) " (6 " i6) 

In view of (2.19) this implies the determinantal equation (cf., after squaring, 
Forrester 2006, Eq. (1.22)) 

D 1 (z; (s r oo))^ = \ (det (l - z{2 - z)V M \ h2[s ^ (l + s J^-^ 

+ det(l+ y Z (2-z)V Ai f L2(S/Co) ) (l-y^)) (6-27) 

for all s 6 R and z in the complex domain of analyticity that contains Z — 1. 
This formula thus paves, following the arguments of the proof of Theorem 6.1, 
an elementary road to the Painleve representation (Dieng 2005, Eq. (1.2.1)) of 
Di(z; (s, 00)) (see also Forrester 2006, Eq. (1.17)). 

Based on the recursion formulae (6.17) and the numerical methods of Section 4 
we calculated the distribution functions Fi(k;s) of the A:-th largest level in the 
soft edge scaling limit of GOE — note that because of (6.14) there is no need 
for a separate calculation of the corresponding distributions F^{k;s) for GSE. 
The corresponding densities are shown, for k = 1, ...,6, in Figure 6. a. These 
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calculations are accurate to the imposed absolute tolerance of 5 ■ 10~ 15 . Taking 
our numerical solutions as reference, Figure 6.b shows the absolute error of the 
Painleve II based numerical evaluations by Dieng (2005), which are available for 
k — 1, ... ,4. Perfectly visible are the points «+ = —10 and b- = 6 where Dieng 
chose to switch from an asymptotic formula to the BVP solution of Painleve II, see 
Section 3.2.1. 

7. The k-th Smallest Level at the Hard Edge: LOE and LSE 

In this section we derive determinantal formulae for the cumulative distribution 
functions Fp A (k;s) of the fc-th smallest level in the hard edge scaling limit of LOE, 
LUE, and LSE with parameter a. It turns out that these formulae have the same 
algebraic structure as those developed for the A:-th largest level in the soft edge 
scaling limit. Therefore, though more concise, we proceed step-by-step in parallel 
to the arguments of the preceding section. 

By the underlying combinatorial structure we have 

k-1 

f pa i k ' s ) = 1 - E %« (/' s )' (7-1) 

;=0 

where we briefly write 

Ep A (k;s) =£{ 5 hard) (/c;(0,s), a ). (7.2) 
The integral representation 

M*,y) = \ [ J* ( Vfr) J* ( % (7.3) 

of the Bessel kernel implies, by introducing the kernels 

T SA {x,y) = —Ja (Vsxy) , V K (x,y) = -/« (y/xy) , (7.4) 
the factorization 

, 2 



det [I - zK a r L2(0/S) J = det ^ - z [T SA \ l2{0i1) 

= det (l - \/zT SA tx,2( 0/1 )J ■ det (l + \/zT SA \ L 2 {0/1) 



= det [1 - Vz V a r L 2(o,vV) J ■ det {! + Vz V« \ L ^0,^)J (7-5) 

that is valid for the complex cut plane z £ C \ (— 00, 0]. Now, upon introducing the 
functions 



(7-6) 

2=1 



we obtain from (2.39) 



k 

E2 A (k;s) = J^E +A (j;s)E- A (k - j;s). {7.7) 

j=o 
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Desrosiers and Forrester (2006, Prop. 1) proved that 

E llt -i(0;s) = E +< «(0;s) (7.8a) 

E 4 , a+ i(0;s) = i(£ + , a (0;s)+£_, ff (0;s)). (7.8b) 
The last generalizes to 

E iA+1 {k;s) = \{E +A {k;s)+E„ A {k;s)) (k = 0, 1,2, . . .), (7-9) 
or, equivalently, to the generating function 

00 

J2E iA+1 (k;s)(l-z) k 

k=0 

= i (det(l-VzV»r t2(0V5) ) +det(7 + V5VJ L 2 (0 ^ ) )) (7.10) 

that holds for all s G (0, 00) and z in the complex domain of analyticity that 
contains 2 = 1. A proof of (7.10), using Painleve representations, was given by 
Forrester (2006, Eq. (1.38)). 

It remains to discuss the LOE. To this end we use the interrelationships between 
LOE, LUE, and LSE found by Forrester and Rains (2001, Thms. 4.3/5.1), 15 which 
can symbolically be written in the form 

LSE nA+ i = even (hOE 2n+l , (7- lla ) 

LUE„ 

,a — even ( LOE^ «_i U LOE^^ «-i j . (7.11b) 

Here, we write LOE,, A for the n-dimensional LOE with parameter a, etc. Otherwise, 
these symbolic equations have the same meaning as (6.13a) and (6.13b). For the 
same purely combinatorial reasons as in the preceding section we hence get 

F^+i(*;s)=.F 1#e -i(2*;s) (k = 1,2,3, .. .), (7.12) 

or, equivalently by (7.1) and (7.9), 

E 1/2 -i(2fc;s) + E 1!£ -i(2fc + l;s) = \{E +A {k;s) + E- A (k;s)), (7.13) 

as well as (see also Forrester 2006, Cor. 4) 

2k 2k+l 

= J^E lX -i(j;s)E lx -i(2k-j;s)+ E i 2=1 0'' s ) E i 2=1 ( 2k + 1 ~^ s )- (?-H) 

7=0 ' 2 ' 2 j=o ' 2 ' 2 

Since Theorem 6.2 was in fact just addressing the solution of a specific system 
of functional equations, we obtain the same result here because the system (7.7), 
(7.8a), (7.13), and (7.14) possesses exactly the algebraic structure considered there. 



I ^Again, that is why we have chosen the weight functions scaled as in Forrester and Rains (2001). 
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That is, we get the solution 

( 2 h 



' '_(?) 
7=0 



E 1/!S -i(2fc;s) = E+, a (fc;s) - £ 2 2;+i ^ + ( 2A: ~ 2 / ~ ^ s )' (7-i5») 



E 1#fi _x(2jfc + l; S ) = g±^M+jz^M - £l ^( 2 fc;s), ( 7 . 15 b) 
and, completely parallel to (6.27), the corresponding generating function 



+ det^+ A /z(2-z)y a r i2(0 , vi) ) (7 .i6) 

for all s 6 (0, 00) and 2 in the complex domain of analyticity that contains z = 1. 
Remark 7.1. Given the large order asymptotics (Olver 1974, Eq. (9.5.01)) 

/,(« + ^ 1/3 ) = 2 1/3 a - 1/3 Ai(-2 1/3 + 0(0 (« -> 00), (7.17) 
which holds uniformly for bounded £, we get, by changing variables subject to 



x = ^u 2 - 2tt(a/2) 1 / 3 £, y = J a 2 - 2a(a/2) 1 / 3 ^, (7.18) 
the kernel approximation 

V a (x,y)^ = -V Ai (Z, n )+0(a~ 2/3 ) (a -> 00), (7.19) 

uniformly for bounded £, 77. Since y = is mapped to r/ = (a/2) 2 / 3 — > 00 (a — > 00) 
we therefore obtain, observing the fast decay of the Bessel and Airy functions, the 
operator approximation (in trace class norm) 

^11(0, ^-2^/2)^ )=^^) +0(«" 2/3 ) (« 00). (7.20) 

This approximation implies the limit of the associated Fredholm determinants, 



Jim det - zV K r L2 ( oV , 2 _ 2|t(<t/2)1/3s ) ) = det [I - zV M y {sfia) ) , (7.21) 
or, equivalently 

Jim E± A (k;a 2 -2a{a/2) 1/3 s) =E±(k;s) (k = 0, 1,2, . . .). (7.22) 

Plugging this limit into (7.7), (7.9), and (7.15) yields, by (6.7), (6.9), and (6.17), the 
hard-to-soft transition (see also Borodin and Forrester 2003, §4) 

Jim E MW (fc ;a 2 -2 a (a/2) 1/3 s) = E^(fc;s) (7.23) 

with 

|>-l)/2, = 1, 

<P) = y, p = 2, (7-24) 

[tx + 1, 0=4. 
This transition is a further convenient mean to validate our numerical methods. 
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8. Matrix Kernels and Examples of Joint Probability Distributions 

8.1. Matrix Kernels. Bornemann (2010a, §8.1) showed that the quadrature based 
approach to the numerical approximation of Fredholm determinants that we 
described in Section 4.1, can fairly easily be extended to matrix kernel determinants 
of the form 



(K u ■■■ K 



l(z) =det \l-z 



IN 



^ 2 (/i)ffi-ffiL 2 (/ N ) ' ( 8 ' 1 ) 



\K N1 ■ ■ ■ K NN/ 

where }\, ■ ■ ■ ,Jn are open intervals and the smooth matrix kernel generates a trace 
class operator on L 2 (/i) ® ■ ■ ■ ® L 2 (/n). By taking an m -point quadrature rule of 
order m with nodes Xu 6 /; and positive weights Wu, written in the form 

m n 

£ WjjfiXij) « / f{x) dx (t = 1 N), (8.2) 

j=l J h 

we approximate d(z) by the N ■ m -dimensional determinant 

/ (A n ■■■ A m \\ 
d m (z) =det I-z : : (8.3a) 

V \A N1 ■■■ A NN ) ) 
with block entries that are the m x m matrices given by 

( A ij)p,q = ™]p 2 Kij(xi P , x jq) w ]q 2 (p = l,...,m;q = l,...,m). (8.3b) 

If the kernel functions Kjj(x,y) are analytic in a complex neighborhood of /, x ]p 
one can prove exponential convergence (Bornemann 2010a, Thm. 8.1): there is a 
constant p > 1 such that 

d m (z) - d{z) = 0(p- m ) (m ->■ 00), (8.4) 

locally uniform in 2 £ C. The results of Section 4.2-4.5 apply then verbatim. 
This approach was used to evaluate the matrix kernel determinant (2.16) for the 
numerical checks of the fundamental equation (6.10) in Section 6. 

8.1.1. An example: the joint probability distribution ofGUE matrix diffusion. Prahofer 
and Spohn (2002) proved that the joint probability of the maximum eigenvalue of 
GUE matrix diffusion at two different times is given, in the soft edge scaling limit, 
by the operator determinant 

F(A 2 (t) < s t ,A 2 (0) ^ s 2 ) = det ^ - r L 2 (siHfflL 2 (s2/00) J (8.5a) 

with kernel 

(rOQ 
J o e -£ f Ai(x + £)A% + £)d£ (t>0), 
lo ( 8 -5 b ) 
e ^ f Ai(x + f)Ai(y + 0^ (*<0). 



(Note that Kq = Kj^, see (6.3).) This expression is directly amenable to our numer- 
ical methods; Figure 7.a shows the covariance function that was calculated this 
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a. cov(yt 2 (f)/ -42(0)) for GUE matrix diffusion b. CDF of the sum of largest two levels in GUE 

Figure 7. Left: plot (blue line) of the covariance cov(^4 2 (0' .4 2 (0)) 
of the maximum eigenvalue of GUE matrix diffusion at two dif- 
ferent times (soft edge scaling limit), calculated from the joint 
probability distribution (8.5) to 10 digits accuracy (absolute toler- 
ance 5 ■ 10~ n ). The red dots show the data obtained from a Monte 
Carlo simulation with matrix dimensions m — 128 and m — 256 
(there is no difference visible between the two dimensions on 
the level of plotting accuracy). The dashed, green lines show the 
asymptotic expansions (8.6) and (8.7). Right: plot (blue line) of the 
CDF (8.24) of the sum of the largest two levels of GUE (soft edge 
scaling limit), calculated from the joint probability distribution 
(8.22) to 10 digits accuracy (absolute tolerance 5 ■ 10 -11 ). To com- 
pare with, we also show (red, dashed line) the convolution (8.25) 
of the corresponding individual CDFs from Figure 4-b; because of 
statistical dependence, there is a clearly visible difference. 



way. The results were cross-checked with a Monte Carlo simulation (red dots), and 
with the help of the following asymptotic expansions (dashed lines): for small t 
with the expansion (Prahofer and Spohn 2002, Hagg 2008) 

cov{A 2 {t),A 2 {Q)) =var(F 2 ) -f + 0(f 2 ) (t ->• 0), (8.6) 

where F 2 denotes the Tracy-Widom distribution for GUE (the numerical value of 
var(F 2 ) can be found in Table 4); for large t with the expansion (Widom 2004, Adler 
and van Moerbeke 2005) 

cav{A 2 {t),A 2 [Q)) = r 2 + cf" 4 + 0(r 6 ) (f 00), (8.7) 

where the constant c = —3.542 ■ ■ • can explicitly be expressed in terms of the 
Hastings-McLeod solution (1.8) of Painleve II. 

Similar calculations related to GOE matrix diffusion, and their impact on "ex- 
perimentally disproving" a conjectured determinantal formula, are discussed in a 
recent paper by Bornemann, Ferrari and Prahofer (2008). 

Remark. In a masterful analytic study, Adler and van Moerbeke (2005) proved 
that the function G(t,x,y) = logP(^2(0 ^ x,<A.2(0) ^ J/) satisfies the following 
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nonlinear 3rd order partial differential equation (together with certain asymptotic 
boundary conditions): 



t. 



dt \ dx 2 



dy 2 



d 3 G ( d 2 G d 2 G d 2 G 



dx 2 dy \ dy 2 dxdy dx 2 



y-t' 



d 3 G fd 2 G _9 2 £_9 2 G_ X+ _ t 2 
dy 2 dx \ dx 2 dxdy dy 2 X 



+ 



fd 3 G d _ d 3 G d \ f d d 
\ dx 3 dy dy 3 dx) \dx dy 



G. (8.8) 



The reader should contemplate a numerical calculation of the covariance function 
based on this PDE, rather than directly treating the Fredholm determinant (8.5) as 
suggested in this paper. 

8.2. Operators Acting on Unions of Intervals. Determinants of integral operators 
K acting on the space L 2 (Ji U ■ ■ ■ U /n) of functions defined on a union of mutually 
disjoint open intervals can be dealt with by transforming them to a matrix kernel 
determinant (Gohberg, Goldberg and Krupnik 2000, Thm. VI.6.1): 



det(J- Z Kf L2(/iU ... u/N) 



det 



V 



K 



IN 



\ 



li 2 (/l)ffi-®I. 2 (7N) 



(8.9) 



where JC« : L 2 (Jj) — > L 2 (/,) denotes the integral operator induced by the given 
kernel function K(x,y), that is, 



Kiju(x) = / K(x,y)u(y)dy (x £ /,). 
JJj 



(8.10) 



More general, along the same lines, with X] denoting the characteristic function of 
an interval / and Zj E C, we have 



N 



det [I- [KjTzjXL )\ L 2 {J 



det 



V 





i. 2 (/i)e-®L2(/ N ) 



(8.11) 



where the : L 2 (Jj) — > L 2 (Ji) denote, once more, the integral operators defined 
in (8.10). To indicate the fact that the operators Kjj share one and the same kernel 
function K(x, y) we also write, by "abus d'langage", 







(K ■ 




det 


h 




1 











(8.12) 
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for (8.9) and 

/ z-[K ■ ■ ■ zpjK^ 



det 



I 



V 



\ztK ■ ■ ■ ZfjKj 

for (8.11). Clearly, in both of the cases (8.9) and (8.11), we then apply the method 
of Section 8.1 for the numerical evaluation of the equivalent matrix kernel deter- 
minant. 

8.3. Generalized Spacing Functions. The determinantal formulae (2.11), (2.14), 
(2.15), (2.16), (2.36), and (2.39) have all the common form 

f_ l)k d k / \ 
E{k; }) = P (exactly k levels lie in /) = det (l -zK\ L2{;) ) , (8.14) 

z=l 

which is based on an underlying combinatorial structure that can be extended 
to describe, for a multi-index a 6 and mutually disjoint open intervals h 
(j — 1,..., N), the generalized spacing function 

E(a; ]\, . . . , Jn) = P(exactly ctj levels lie in Jj, j = 1, . . . , N) (8.15a) 

by the determinantal formula (see Tracy and Widom (1993a), Thm. 6, or (1998), 
Eq.( 4 .i)) 



EO;Ji,...,/n) 



N 



. (8.15b) 

By the results of Section 8.2 and 4.3, such an expression is, in principle at least, 
amenable to the numerical methods of this paper. However, differentiation with 
respect to / different variables Zj means to evaluate an Z-dimensional Cauchy inte- 
gral of a function that is calculated from approximating a Fredholm determinant. 
This can quickly become very expensive, indeed. 

8.3.1. Efficient numerical evaluation of the finite-dimensional determinants. Let us 
briefly discuss the way one would adapt the methods of Section 4.2 to the current 
situation. We will confine ourselves to the important case of a multi-index a which 
has the form a = (k, 0, . . .,0). If we are to compute the derivatives by a Cauchy 
integral, we have to evaluate a determinant of the form 16 

d(z) = det (I - (zA \ B)) A G R m 'P,B £ R m ' m ~P (8.16) 

for several different arguments zeC. By putting 

T = I-(0\B), A = T _1 (A 1 0), (8.17) 

we have I — (zA \ B) = T(I — zA) and hence the factorization 

d(z) = det(T) • det(J - zA) (8.18) 

to which the results of Section 4.2 apply verbatim. 

l6 Which is understood to be an approximation of the determinant in (8.15), expressed as the matrix 
kernel determinant (8.11). 
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8.3.2. An example: the joint distribution of the largest two eigenvalues in GUE. Let us 
denote by X\ the largest and by A 2 the second largest level of GUE in the soft edge 
scaling limit (2.4). We want to evaluate their joint probability distribution function 

F{x,y) = P(A a ^x,A 2 (8.19) 

Since Ai ^ x ^ y certainly implies A2 ^ y, we have 

F(x,y) = P(Ax ^ x) = F 2 {x) (x < y). (8.20) 

On the other hand, if x > y, the open intervals (y, x) and (x, 00) are disjoint and 
we obtain for simple combinatorial reasons 

F(x,y) = E < soft) (0, 0; (y,x), (x, 00) ) + E ( soft) (1, 0; (y, at), (x, 00) ) 

= £( soft) (0; (y,oo)) + Ef ft) (l,0; (y,x), (*,oo)) 

= F 2 (y) + E 2 ;soft) (l,0;(y,x),(x,co)) (x > y). (8.21) 
By (8.15) and (8.11) we finally get 

F(x,y) 

'F 2 {x) (x^y), 




f >"/> - ^ det ( I - [ ^ ^ j rL2 (y ,, )eL2(X/0o) 



(x >y). 



(8.22) 



/ 2=1 

We have used this formula to calculate the correlation of the largest two levels as 

p(A 1 ,A 2 ) = 0.50564 723159...; (8.23) 

where all the 11 digits shown are estimated to be correct and the run time was 
about 16 hours using hardware arithmetic. 17 So, as certainly was to be expected, 
the largest two levels are statistically very much dependent. Figure y.b plots the 
CDF of the sum of the largest two levels, 

/oo 
3iF(x,s-x)dx, (8.24) 
-CO 

and compares the result with the convolution of the individual CDFs of these 
levels, that is with the s-dependent function 

/oo 
F 2 , (l;x)F 2 (2;s-x)dx. (8.25) 
-00 

Fhe clearly visible difference between those two functions is a further corollary of 
the statistical dependence of the largest two level. 

Yet another calculation of joint probabilities in the spectrum of GUE (namely, 
related to the statistical independence of the extreme eigenvalues) can be found in 
Bornemann (2010b). 



I7 So, this example stretches our numerical methods pretty much to the edge of what is possible. 
Note, however, that these numerical results are completely out of the reach of a representation by 
partial differential equations. 
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Table 5. Toolbox commands for kernel functions K(x,y). If the 
value of K(x,y) is defined as an integral, there is an additional ar- 
gument m to the command that assigns the number of quadrature 
points to be used. 



kernel 


formula 


command 


vectorization mode 


Y (v m\ 


(1.2) 


sine (pi*(x-y)) 


' grid ' 


V&Ax v) 

v Al \^-f y } 


(6 4.) 


airv( (x+v) /2) 11 


'grid' 


V a (x,y) 


(74) 


besselj (alpha, sqrt (x. *y) ) /2 


'grid 


K„(x,y) 


(2.12) 


Hermit eKernel (n , x , y ) 


'outer' 


K M (x,y) 


(1.6) 


AiryKernel (x , y ) 


'outer' 


K nA (x,y) 


(2-37) 


LaguerreKernel (n , alpha , x , y ) 


'outer' 


K a (x,y) 


(2.39c) 


BesselKernel (alpha, x,y) 


'outer' 


S(x,y) 


(2.17a) 


F4MatrixKernel(x,y,m, 'SN') 


'outer' 


S*(x,y) 


(2.17a) 


F4MatrixKernel(x,y,m, 'ST') 


'outer' 


SD(x,y) 


(2.17b) 


F4MatrixKernel(x,y,m, 'SD') 


'outer' 


IS(x,y) 


(2.17c) 


F4MatrixKernel(x,y,m, 'IS') 


'outer' 


Kt(x,y) 


(8.5b) 


Airy2ProcessKernel (t , x , y , m) 


'outer' 



9. Software 

We have coded the numerical methods of this paper in a Matlab toolbox. (For 
the time being, it can be obtained from the author upon request by e-mail. At a 
later stage it will be made freely available at the web.) In this section we explain 
the design and use of the toolbox. 

9.1. Low Level Commands. 

9.1.1. Quadrature rule. The command 
» [w,x] = ClenshawCurtis(a,b,m) 

calculates the ffz-point Clenshaw-Curtis quadrature rules (suitably transformed if 
a = —00 or b — 00). The result is a row vector w of weights and a column vector x 
of nodes. This way, the application (4.1) of the quadrature rule to a (vectorizing) 
function f goes simply by the following command: 

» w*f(x) 

Once called for a specific number m, the (untransformed) weights and nodes are 
cached for later use. As an alternative the toolbox also offers Gauss-Jacobi quadra- 
ture, see Section A.i for its use in the context of algebraic kernel singularities. 

9.1.2. Kernels and vectorization modes. The numerical approximation of Fredholm 
determinants, by (4.3) or (8.3), requires the ability to build, for given m-dimensional 
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column vectors x and y, the m x m matrix 

A = (K(x„y ; -))"; = i- 

For reasons of efficiency we make a meticulous use of Matlab's vectorization 
capabilities. Depending on the specific structure of the coding of the kernel 
function, we distinguish between two vectorization modes: 

(1) 'grid': The matrix A is built from the vectors x, y and the kernel function 
K by the commands 

» [X,Y] = ndgrid(x,y) ; 
» A = K(X,Y) ; 

(2) ' outer ' : The matrix A is built from the vectors x, y and the kernel function 
K by the commands 

» A = K(x,y); 

Table 5 gives the commands for all the kernels used in this paper together with 
their vectorization modes. 

9.1.3. Approximation ofFredholm determinants. Having built the matrix A the ap- 
proximation (4.3) is finally evaluated by the following commands: 

» w2 = sqrt (w) ; 

» det(eye(size(A))-z*(w2'*w2) . *A) 

9.1.4. Example. Let us evaluate the values Fj(0) and F2(0) of the Tracy-Widom 
distributions for GOE and GUE by these low level commands using m — 64 
quadrature points. The reader should observe the different vectorization modes: 

» m = 64; [w,x] = ClenshawCurtis (0 , inf ,m) ; w2 = sqrt(w); 

» [xi,xj] = ndgrid(x,x); 

» Kl = @(x,y) airy((x+y)/2)/2; 

» F10 = det(eye(m)-(w2'*w2) . *K1 (xi ,xj ) ) 

F10 = 0.831908066202953 

» KAi = OAiryKernel; 

» F20 = det(eye(m)-(w2'*w2) .*KAi(x,x)) 

F20 = 0.969372828355262 

A look into Prahofer's (2003) tables teaches that both results are correct to one unit 
of the last decimal place. 

9.2. Medium Level Commands. The number of quadrature points can be hidden 
from the user by means of the automatic error control of Section 4.4. This means, 
we start thinking in terms of the limit of the approximation sequence, that is, in 
terms of the corresponding integral operators. This way, we evaluate the operator 
determinant 

det(I-zK\ L2{J) ) 
for a given kernel function K(x,y) by the following commands: 
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» K.ker = (3(x,y) K(x,y); k.mode = vectorizationmode ; 

» Kop = op(K, J) ; 

» [val, err] = detlm(Kop,z) ; 

(The argument z may be omitted in detlm if z = 1.) Here, val gives the value of 
the operator determinant and err is a conservative error estimate. The code tries 
to observe a given absolute tolerance tol that can be set by 

» pref ('tol' ,tol) ; 

The default is 5 ■ 10 -15 , that is, tol = 5e-15. The results can be nicely printed in 
a way such that, within the given error, either just the correctly rounded decimal 
places are displayed (printmode = 'round'), or the correctly truncated places 
(printmode = 'trunc'): 

» pref ( 'printmode ' , printmode) ; 
» PrintCorrectDigits(val,err) ; 

For an integral operator K, the expressions 



(-1)* d k 
k\ dz 1 



det{I- zK) 



2=1 



(-1)* d k 
k\ dz 1 



det(I - VzK) 



(-1)* 



2=1 



det(J - zK) 



k\ dz k 

are then evaluated, with error estimate, by the following commands: 

» [val, err] = dzdet(K.k); 

» [val, err] = dzdet(K,k,@sqrt) ; 

» [val, err] = dzsqrtdet (K,k) ; 



2=1 



9.2.1. Example 9.1.4 revisited. Let us now evaluate the values F\(0) and £2(0) of the 
Tracy- Widom distributions for GOE and GUE by these medium level commands. 

» pref ( 'printmode ', 'trunc ') ; 

» Kl.ker = @(x,y) airy ( (x+y) /2) /2 ; Kl.mode = 'grid'; 
» [val, err] = detlm(op(Kl , [0, inf ] ) ) ; 
» PrintCorrectDigits (val , err) ; 

0.83190806620295, 

» KAi.ker = SAiryKernel; KAi.mode = 'outer'; 
» [val, err] = detlm(op(KAi , [0 , inf ] ) ) ; 
» PrintCorrectDigits (val, err) ; 

0.96937282835526_ 

So, the automatic error control supposes 14 digits to be correct in both cases; a 
look into Prahofer's (2003) tables teaches us that this is true, indeed. 
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9.2.2. Example: an instance of equation (6.10). We now give the line of commands 
that we used to experimentally check the truth of the determinantal equation (6.10) 
before we worked out the proof. For a single instance of a real value of s and a 
complex value of z we obtain: 

» s = -1.23456789; z = -3.1415926535 + 2 . 7182818284i ; 

» 

» Kll.ker = @(x,y,m) F4MatrixKernel (x,y ,m, 'SN' ) ; Kll.mode = 'outer'; 
» K12.ker = @(x,y,m) F4MatrixKernel(x,y ,m, 'SD ' ) ; K12.mode = 'outer'; 
» K21.ker = @(x,y,m) F4Matr ixKernel (x,y,m, 'IS') ; K21.mode = 'outer'; 
» K22.ker = @(x,y,m) F4Matr ixKernel (x , y , m , ' ST ' ) ; K22.mode = 'outer'; 
» K = op({Kll K12; K21 K22> , { [s , inf ] , [s , inf ] }) ; 
» vail = sqrt(detlm(K,z)) 

vail = 1.08629916321436 - . 0746712169305511i 

» Kl.ker = @(x,y) airy ( (x+y) /2) /2 ; Kl.mode = 'grid'; 
» K = op(Kl, [s.inf]) ; 

» val2 = (detlm(K,sqrt(z)) + detlm(K, -sqrt (z) ) ) /2 

val2 = 1.08629916321436 - . 0746712169305508i 
» dev = abs(vall-val2) 

dev = 5.23691153334427e-016 

This deviation is below the default tolerance 5 • 10~ 15 which was used for the 
calculation. 

9.3. High Level Commands. Using the low and medium level commands we 
straightforwardly coded all the functions that we have discussed in this paper 
so far. Table 6 lists the corresponding commands. The reader is encouraged to 
look into the actual code of these commands to see how closely we followed the 
determinantal formulae of this paper. 

9.3.1. Example 9.2.1 revisited. Let us evaluate, for the last time in this paper, the 
values Fi(0) and F2(0) of the Tracy-Widom distributions for GOE and GUE, now 
using those high level commands. 

» pref ( 'printmode ' , ' trunc ' ) ; 

» [val, err] = F(l ,0) ; 

» PrintCorrectDigits (val , err) ; 

0.83190806620295. 

» [val, err] = F(2,0) ; 

» PrintCorrectDigits (val, err) ; 

0.96937282835526_ 

9.3.2. Example: checking the k-level spacing functions against a constraint. An appropri- 
ate way of checking the quality of the automatic error control goes by evaluating 
certain constraints such as the mass and mean given in (2.8): 
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Table 6. Toolbox commands for all the probability distributions 
of this paper. A call by val = E( . . . ) etc. gives the value; with 
[val , err] = E ( . . . ) etc. we get the value and a conservative 
error estimate. 



f n n r"H nn 


Hpfinino" formnlpip 

\_1 1 1 1 Ll± L i. Ul 11L L4.10.tr 


rnmmanH 

Lull Lll LQ1 l\_l 


infprval J — 1 Si ) 

J.± l Lt 1 V C\ 1 1 1 j 1 , iZ> J J 

interval / = (s, 00) 
interval / = (— 00, s) 




T = r=i =91 
J = [s.inf] 
J = [-inf,s] 


F^((k nv h M 
E< bulk) (fc;J) 
^((fcO);/^) 
4 soft) (fc;7) 

2_ \ ' 3 ; 

4 oft \(k,0);h,} 2 ) 
F(x,y) 


(2.1l) 

(2.14) 
(8.15) 
(2.15) 
(8.15) 
(8.22) 


E(2,k, J,n) 

FC9 ric m -T T1 T9> Til 

E(2,k, J, 'bulk') 

E(2, [k,0] ,4JJ1 ) J2}, 'bulk') 

E(2,k,J, 'soft') 

E(2, [k,0] ^Jl,.^}, 'soft') 

F2 Joint (x,y) 


Ei S ° tt] (k;J) 


(2.16) 


E(4,k, J, 'soft' , 'MatrixKernel') 


E^(k;J,K) 
EiiL((k,Q);h h.ct] 
Ef mA \k;],a) 
E^\{k,Q);h,h,oc) 


(2.36) 
(8 is) 

(2.39) 

(8.15) 


ECLUE' ,k, J.n.alpha) 

E('LUE' , [k,0] ,{J1, J2},n, alpha) 

E(2,k, J, 'hard' .alpha) 

E(2, [k,0] ,{J1,J2}, 'hard' .alpha) 


E+(k;s) 
E-(hs) 
E/j(fc;s) 


(5-6) 
(5-6) 

(2-5), (5-7), (5-8), (5-9) 


E('+' ,k,s) 
E('-' ,k,s) 
E(beta,k,s) 


£+M 
£_(fc;s) 


(6.6) 
(6.6) 

(6.2), (6.7), (6.9), (6.17) 

(2.9), (6.1) 

(2.10) 


E('+' ,k,s, 'soft') 
E('-',k,s, 'soft') 
E(beta,k,s, 'soft') 
F(beta,k,s) 
F(beta,s) 


E + , a (fc;s) 
E-, a (fc;s) 
Ep A (k;s) 
Ef, A {k;s) 


(7-6) 
(7-6) 

(7-2), (7-7), (7-9)- (7-15) 
(7-i) 


E('+' ,k,s, 'hard' .alpha) 
E('-',k,s, 'hard' .alpha) 
E (beta , k , s , ' hard ' , alpha) 
F(beta,k,s, alpha) 
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» pref ( 'printmode ', 'round' ) ; 
» s = 2.13; beta = 1; 

» mass = 0; errmass = 0; mean = 0; errmean = 0; 

» M = 10; for k=0:M 

» [val.err] = E(beta,k,s); 

» mass = mass+val; errmass = errmass+err; 

» mean = mean+k*val; errmean = errmean+k*err ; 

» end 

» PrintCorrectDigits (mass , errmass) ; 



1.0000000000000. 



» PrintCorrectDigits (mean, errmean) ; 
2. 130000000000__ 

The results of (2.8) are perfectly matched. The reader is invited to repeat this 
experiment with a larger truncation index for the series. 

9.3.3. Example: more general constraints. The preceding example can be extended to 
more general probabilities E(k; J) that are given by a determinantal expression of 
the form (8.14), that is, 



(9-i) 



z=l 



E(fc;/)=^ d? det(l-zKr L2(/) ) 
for some trace class operator K\ L 2gy Expanding the entire function 

d(z) =det(j-zJCr L 2 (/) ) (9-2) 



into a power series at z = 1 yields 

E = E ^^(1) = d(o) = 1, (9.3a) 



A; 1 

k=0 k=0 



£*e(* ; j) = - E Sr Ld(fc+1) ( 1 ) = ~ d 'W = tr ; (9.3b) 

k=0 k=o K - 



both of which have a probabilistic interpretation (see Deift 1999b, p. 119). Now, for 
the Airy kernel we get 



tr (^AiL2( s ,oo)) = J ^AiO,*) dx 



— 3 

with the specific value (for s = 0) 



i(2s 2 Ai(s) 2 -2sAi'(s) 2 - Ai(s)Ai'(s)), (9.4) 



te (^ fe ^)) = 9niki)' {9 ' 5) 

Now, let us check the quality of the numerical evaluation of (9.3) (and the automatic 
error control) using this value. 
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» pref ( 'printmode ', 'round' ) ; 
» s = 0; beta = 2; 

» mass = 0; errmass = 0; mean = 0; errmean = 0; 
» M = 3; for k=0:M 

» [val, err] = E(beta,k,s, 'soft ' ) ; 

» mass = mass+val; errmass = errmass+err; 

» mean = mean+k*val; errmean = errmean+k*err ; 

» end 

» PrintCorrectDigits (mass , errmass) ; 

1.000000000000__ 
» PrintCorrectDigits (mean, errmean) ; 

0.030629383079___ 
» l/9/gamma(l/3)/gamma(2/3) 

0.0306293830789884 

The results are in perfect match with (9.3): they are correctly rounded, indeed. The 
reader is invited to play with the truncation index of the series. 

9.3.4. Example: calculating quantiles. Quantiles are easy to compute; here come the 
5% and 95% quantiles of the Tracy-Widom distribution for GOE: 

» Flinv = vec(@(p) (fzero(@(s) F(l , s) -p, 0) ) ) ; 
» Flinv([0.05 0.95]) 

-3 . 18037997693773 . 979316053469556 

9.4. Densities and Moments. Probability densities and moments are computed 
by barycentric interpolation in m Chebyshev points as described in Section 4.5. 
This is most conveniently done by installing the functionality of the chebfun 
package (Trefethen, Pachon, Platte and Driscoll 2008). Our basic command is then, 
for a given cumulative distribution function F(s): 

» [val , err , supp , PDF , CDF] = moments(@(s) F(s),m); 

If one skips the argument m, the number of points will be chosen automatically. The 
results are: val gives the first four moments, that is, mean, variance, skewness, and 
kurtosis of the distribution; err gives the absolute errors of each of those moments; 
supp gives the numerical support of the density; PDF gives the interpolant (4.13) 
of the probability density function F'(s) in form of a chebfun object; CDF gives the 
same for the function F (s) itself. If one sets 

» pref ( 'plot true) ; 

a call of the command moments will plot the PDF and the CDF in passing. 

9.4.1. Example: the first four moments of the Tracy-Widom distributions for GOE, GUE, 
and GSE. The following fills in the missing high-precision digits for the GSE in 
Tracy and Widom (2008, Table 1); Figure 8 is automatically generated in passing: 
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Figure 8. Probability density functions Fg(s) (PDF, blue) and 
cumulative distribution functions (s ) (CDF, green) of the Tracy- 
Widom distributions (2.10) for GOE (j6 = 1), GUE (/S = 2), and 
GSE (/3 = 4); larger f> go to the left. Fhis plot was automatically 
generated by the commands in Example 9.4.1. Compare with 
Tracy and Widom (1996, Fig. 1) and Dieng (2005, Fig. 1.1). 

» pref ( 'printmode ' , ' trunc ' ) ; pref ( 'plot ' ,true) ; 
» for beta =[12 4] 

» [val.err] = moments (@(s) F(beta,s), 128); 
» PrintCorrectDigits (val , err) ; 
» end 

-1.2065335745820. 1 . 607781034581__ 0.29346452408 0.1652429384. 

-1.771086807411__ . 8131947928329__ . 224084203610___ 0.0934480876. 
-2.306884893241__ . 5177237207726__ 0.16550949435 0.0491951565. 



9.4.2. Example: checking the k-level spacing functions against a constraint. In the final 
example of this paper we check our automatic error control in dealing with 
densities. We take the level spacing densities defined in (2.6) and evaluate the 
integral constraints (2.7): 

» pref ( 'printmode ', 'round' ) ; beta = 1; M = 10; [dom.s] = domain (0,M); 

» E0 = chebfun(vec(@(s) E(beta, , s) ) ,dom) ; 

» El = chebfun(vec(@(s) E(beta, 1 , s) ) ,dom) ; 

» E2 = chebfun(vec(@(s) E(beta, 2 , s) ) ,dom) ; 

» p2 = diff (3*E04-2*E1+E2,2) ; 

» mass = sum(p2) ; errmass = cheberr (p2) ; 

» mean = sum(s.*p2); errmean = errmass*sum(s) ; 

» PrintCorrectDigits (mass , errmass) ; 

1.0000000000 

» PrintCorrectDigits (mean, errmean) ; 

3.000000000 
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Once more, the results of (2.7) are perfectly matched. The reader is invited to 
repeat this experiment with a larger truncation point for the integrals. 

A. Appendices 

A.i. Algebraic Kernel Singularities and Gauss-Jacobi Quadrature. If the kernel 
is not analytic in a complex neighborhood of the interval defining the integral 
operator, the straightforward use of the method of Section 4.1, based on Clenshaw- 
Curtis quadrature, does not generally yield exponential convergence. Since the 
software limits the number of quadrature points (default is a maximum of 256) this 
will result in rather inaccurate results (inaccuracies which are detected, however, 
by the automatic error control). This remark applies in particular to the Laguerre 
kernel K„ A (x,y) and the Bessel kernel K a (x,y) which, for non-integer parameter a, 
exhibit algebraic singularities at x — or y = 0, namely 

K nA (x,y) a x a/2 (x->0), K a {x,y) o< x a/2 (x 0), (A.ia) 

K nA (x,y)<xy a/2 (y->0), K a (x,y) oc y^ 2 (y->0). (A.ib) 

A.I.I. Example: level spacing at the hard edge for non-integer parameter. Hence, for 
example, if we want to evaluate the values 

Ef rd) (l;(0,6),±i) = - jUet(l-zK ± i \ L * m )\ _, 

the method of Section 4.1 gives, using Clenshaw-Curtis quadrature, just 1 digit 
for the stronger singularity a = — j and 6 digits for the weaker singularity a = \ : 

» pref (' quadrature ' ,@ClenshawCurtis) ; pref ( 'printmode ' , 'trunc') ; 
» s = 6; 

» alpha = -0.5; K.ker = @(x,y) BesselKernel (alpha, x, y) ; K.mode = 'outer'; 
» [vall.errl] = dzdet(op(K, [0,s] ) , 1) ; 

» alpha = 0.5; K.ker = @(x,y) BesselKernel(alpha,x,y) ; K.mode = 'outer'; 

» [val2,err2] = dzdet(op(K, [0,s] ) , 1) ; 

» PrintCorrectDigits( [vail val2] , [errl err2] ) ; 

0.8 0.524976 



A. 1.2. Gauss-Jacobi quadrature. This loss of accuracy can be circumvented if we 
switch from Clenshaw-Curtis quadrature to Gauss-Jacobi quadrature (Szego 1975, 
§15.3), which exactly addresses this type of algebraic singularities at the boundary 
points of the interval. Specifically, this quadrature rule addresses functions / on 
the interval [a, b) such that, for oc,fi> —1, the function / that is obtained from 
removing the singularities 

f(x) <x (x - a) a (x — s> a), f(x) oc (b — x)& [x — s> b), (A. 2) 

by 
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is analytic in a complex neighborhood of (a, b). Then, the m-point Gauss-Jacobi 
quadrature provides nodes Xj 6 (a, b) and positive weights Wj such that the approx- 
imation 

m rb 

J^ w jf( x j)~ f(x)dx (A.4) 
j=l Ja 

is exponentially convergent: there is a p > 1 such that the error behaves like 
0(p~ m ) for m — > 00. In our toolbox the Gauss-Jacobi quadrature is called by: 

» [w,x] = GaussJacobi(a,b, alpha, beta, m) 

A. 1.3. Example A.1.1 revisited. The application of the Gauss-Jacobi quadrature to 
the Bessel kernel determinant (2.39) requires the determination of the singular 
exponents that are relevant for the Fredholm determinant. Now, from (A.i) we 
infer that the eigenfunctions 

/ K a {x,y)u K {y)dy = Au A (x) (A.5) 
Jo 

of the Bessel kernel exhibit the same type of algebraic singularity at x — » 0: 

u A (x)otx a/1 (x ->■()). (A.6) 

It turns out that the singularity of K a (x,y)u\(y) at y — >■ governs the behavior of 
the determinant approximation, that is, we should take the singular exponent of 

K a (x,y)u A (y) cxy a (y -> 0). (A.7) 

Indeed, with this choice of the singular exponent we get 14 digits for both cases: 

» alpha = -0.5; pref ( 'quadrature ', @(a,b,m) GaussJacobi (a,b, alpha, 0,m) ) ; 
» K.ker = @(x,y) BesselKernel(alpha,x,y) ; K.mode = 'outer'; 
» [vail, errl] = dzdet(op(K, [0,s] ) , 1) ; 

» alpha = 0.5; pref ( 'quadrature ', @(a,b,m) GaussJacobi (a, b .alpha, 0,m) ) ; 

» K.ker = @(x,y) BesselKernel(alpha,x,y) ; K.mode = 'outer'; 

» [val2,err2] = dzdet(op(K, [0,s] ) , 1) ; 

» PrintCorrectDigits ( [vail val2] , [errl err2] ) ; 

. 86114217058328_ . 52497677921859_ 

All this is most conveniently hidden from the user who can simply call the high 
level command that takes care of choosing the appropriate quadrature rule: 

» alpha = -0.5; [vail, errl] = E(2 , 1 , [0, s] , 'hard' , alpha) ; 
» alpha = 0.5; [val2,err2] = E(2 , 1 , [0, s] , 'hard' , alpha) ; 
» PrintCorrectDigits ( [vail val2] , [errl err2] ) ; 

. 86114217058328_ . 52497677921859_ 

Remark. We have studied this particular example because it provides a cross check 
by the relation 

Ef ard) (fc;(0,s),±i) =E T {k;2V~s/n). (A.8) 
This relation can be obtained from a simple integral transformation that allows us 
to rewrite the determinant of the Bessel kernel for a = ± \ as that of the odd or 
even sine kernel. This way (switching back to Clenshaw-Curtis quadrature which 
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Table 7. Statistical properties of p\{k;s) for various k. Note that 
because of p${k;s) = 2pi(2k + 1,2s) one can read off the values 
for p${k;s) from those for p%(2k + l;s): Just divide the mean by 
two and the variance by four, skewness and kurtosis remain un- 
changed. 



PDF 


mean 


variance 


skewness 


kurtosis 


time 


pi(0;s) 


1 


0.2855306557 


0.6871899889 


0.3712380638 


0.67 


Pi(i; s ) 


2 


0.4163936889 


0.3493968438 


0.0285827332 


1.22 


Pi(2;s) 


3 


0.4974552604 


0.2274144134 


-0.0132956588 


2.59 


pi(3;s) 


4 


0.5556424180 


0.1664568639 


-0.0199468028 


4.24 


Pi(4;s) 


5 


0.6009183521 


0.1304207251 


-0.0200729233 


5.81 


pi(5;s) 


6 


0.63794 46245 


0.10679 47124 


-0.0188407449 


7.43 


pi(6;s) 


7 


0.6692553948 


0.0901832871 


-0.0174319487 


9.46 


Pi(7;s) 


8 


0.6963760657 


0.0779015490 


-0.0161354800 


11.04 


pi(8;s) 


9 


0.7202945046 


0.0684707897 


-0.0150075200 


12.80 


pi(9;s) 


10 


0.74168 65573 


0.06101 25387 


-0.0140407984 


15.26 



is appropriate here) we get values that are in perfect agreement with the ones 
obtained for the Bessel kernel with Gauss-Jacobi quadrature: 

» [vall.errl] = E( '+' , 1 , 2*sqrt (s) /pi) ; 
» [val2,err2] = E( ' - ' , 1 , 2*sqrt (s) /pi) ; 
» PrintCorrectDigits ( [vail val2] , [errl err2] ) 

. 861142170583288 . 524976779218593 

A. 2. Tables of Some Statistical Properties. We provide some tables of statistical 
properties of the A:-level spacing densities pp(k;s), defined in (2.6), and the distri- 
butions Fp(k;s) of the A:-th largest level at the soft edge, defined in (2.9). The tables 
display correctly truncated digits that have passed the automatic error control of 
the software in Section 9. Computing times are in seconds. 

Note that because of the interrelations between GOE and GSE there is no need 
to separately tabulate the values for /3 = 4: first, the interlacing property (6.14) 
gives F^(k;s) = Fi(2k;s). Second, we infer from (5.8) and (5.9) that 

E 4 (k;s) = |E 1 (at-l;2s) + E 1 (2fe2s) + |£i(2fc + l;2s) (A.9) 

which implies 

Pi(k;s) = 2pi(2fc + l;2s). (A.10) 
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Table 8. Statistical properties of p 2 for various k. 



PDF 


mean 


V CI 1 1 CI 1 LV_t- 


skpwnpss 


kurtosis 


time 


P2\v,S) 


i 
i 


fl 1 7QQQ 3877A 




vj.lZ.bby WioU 


U.DO 






f) 94897 77536 


9416743158 


—0 01 494 93984 


1 36 


P2(2;s) 


3 


0.2901698290 


0.15542 00591 


-0.0231740428 


1.44 


p 2 (3;s) 


4 


0.3194435563 


0.1133461773 


-0.0215023114 


1.67 


p 2 (4;s) 


5 


0.3421408054 


0.08871 43069 


-0.0191418388 


2.09 


p 2 (5;s) 


6 


0.3606745961 


0.0726343907 


-0.0171428515 


2.06 


p 2 (6;s) 


7 


0.3763363928 


0.0613508835 


-0.0155525979 


2.19 


P2(7;s) 


8 


0.38989 74631 


0.05301 56552 


-0.01428 79010 


2.41 


p 2 (8;s) 


9 


0.4018551105 


0.04661 73337 


-0.0132681121 


2.56 


p 2 (9;s) 


10 


0.4125486854 


0.04155 73856 


-0.0124320513 


2.74 



Table g. Statistical properties of Fi(fc;s) for various k. Note that 
because of F 4 (A:;s) — F\{7k,s) one can directly read off the values 
for F 4 (fc;s) from those for Fx(2fc;s). 



CDF 


mean 


variance 


skewness 


kurtosis 


time 


Fi(l;s) 


-1.2065335745 


1.6077810345 


0.2934645240 


0.1652429384 


4.59 


Fi(2;s) 


-3.26242 79028 


1.03544 74415 


0.1655094943 


0.0491951565 


12.45 


Fi(3;s) 


-4.8216302757 


0.8223901151 


0.11762 14761 


0.0197746604 


30.04 


Fi(4;s) 


-6.1620399636 


0.7031581054 


0.0923283954 


0.0081606305 


51.24 


h(5;s) 


-7.3701147042 


0.6242523679 


0.0765398210 


0.0024540580 


77.49 


Fi(6;s) 


-8.48621 83723 


0.56700 71487 


0.0656707705 


-0.0007342515 


112.00 




Table io. Statistical properties of F 2 (A:;s) for 


various k. 




CDF 


mean 


variance 


skewness 


kurtosis 


time 


F 2 (l;s) 


-1.7710868074 


0.8131947928 


0.2240842036 


0.0934480876 


1.84 


F 2 (2;s) 


-3.67543 72971 


0.5405450473 


0.1250270941 


0.02173 96385 


5.44 


F 2 (3;s) 


-5.17132 31745 


0.43348 13326 


0.0888080227 


0.0050966000 


10.41 


F 2 (4;s) 


-6.47453 77733 


0.3721308147 


0.0697092726 


-0.0011415160 


17.89 


F 2 (5;s) 


-7.6572422912 


0.33101 06544 


0.0577755438 


-0.0040583706 


25.56 


F 2 (6;s) 


-8.75452 24419 


0.3009494654 


0.04955 14791 


-0.0055998554 


34.72 
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